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IlT  ARSTHACT 


This  p,p«r  presents  the  result!  of  a  study  conducted  to 
develop  algorithms  for  ordering  and  organlilng  data  that  can 
he  presented  In  a  two-dimensional  matrix  form.  The  purpose 
uas  to  osvelop  methods  to  extract  latent  data  patterns,  group¬ 
ings.  and  structural  relationships  which  are  not.  In  general, 
apparent  from  raw  matrix  data. 

Three  distinct  algorithms  were  developed  and  are  present¬ 
ed.  They  are  applied  to  a  variety  of  examples  from  the  social 
and  technical  sciences  which  are  discussed.  The  first  method 
develooed ,  the  Homent  Ordertng  Algorithm,  Is  an  effective  tech¬ 
nique  for  uncovering  and  displaying  a  dominant  univariate  re¬ 
lationship  between  two  sets  of  entitles  that  lie  along  the 
ve-tlcal  and  horiconta!  axes  of  a  matrix.  The  second  method, 
the  Homent  Compression  Algorithm,  Is  designed  to  factor  decom¬ 
posite  matrices  by  proper  reordering.  The  last  method  develop¬ 
ed,  the  Bond  Energy  Algorithm,  was  found  to  be  applicable  to  a 
broader  class  of  problems  than  the  first  two  methods  and  Is  able 
to  efficiently  organwe.  group,  and  Interrelate  data  of  consider¬ 
ably  wore  complex  structure.  The  techniques  developed  are 
applicable  to  c  variety  of  problems  involving  multivariate  deta 
amalysls  and  can  significantly  augment  the  level  of  understand¬ 
ing  and  cOBprehens ion  of  complicated  multivariate  relationships. 
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FOREWORD 


This  paper  presents  the  results  of  a  study  undertaken 
to  develop  methods  for  ordering  a^d  organizing  technical, 
social,  economic  and  other  data  that  can  be  presented  in  array 
form.  The  study  leading  to  the  development  of  this  report 
was  conducted  as  independent  research  at  the  Institute  for 
Defense  Analyses.  The  theory  and  development  of  the 
algorithms  described  in  this  paper  are  the  work  of  members  of 
the  Systems  Evaluation  Division. 
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I.  INTRODUCTION 


Since  the  introduction  of  the  large  digital  computers,  methods  of  multivariate  analysis1 2 
are  being  developed  that  utilize  more  effectively  the  computational  resources  and  character¬ 
istics  of  the  computer  than  some  of  the  more  conventional  and  established  statistical  tech¬ 
niques.  These  new  methods  are  being  employed  because  it  is  now  possible  to  undertake  data 
analysis  problems  in  considerably  greater  detail  than  was  previously  feasible.  A  class  of 
techniques  that  is  able  to  account  for  detailed  individual  relationships  as  well  as  macroscopic 
data  structure  is  exemplified  by  the  cluster-seeking1  methods.  Ball  (Ref.  1)  has  accurately 
pointed  out  that  many  classical  statistical  techniques  depend  heavily  on  statistical  quantities 
estimated  from  the  data  and  that  this  “averaging”  from  the  data  can  sometimes  lead  to 
erroneous  conclusions.  This  is  simply  because  microscopic  variations  in  the  data  cannot,  in 
general,  be  detected  from  the  statistical  quantities  estimated  with  the  result  that  small  but 
significant  information  can  be  overwhelmed  and  even  lost  under  the  pressure  of  larger 
statistical  trends.  Furthermore,  many  of  these  classical  techniques  such  as  principal  component 
analysis  (Ref.  28)  or  factor  analysis  (Ref.  28)  implicitly  assume  data  distributions  that  are  not 
always  present.  Thus,  it  appears  that  there  is  a  definite  need  for  better  direct  analysis 
techniques  so  that  it  is  not  necessary  to  completely  rely  on  functions  of  data  or  on 
assumptions  regarding  their  distribution. 

This  paper  presents  three  new  direct  data  analysis  techniques  that  were  developed  at 
the  Institute  for  Defense  Analyses.  One  of  the  algorithms,  the  Bond  Energy  Algorithm,  shares 
a  few  of  the  same  objectives  as  some  of  the  other  cluster-seeking  techniques  (Refs.  2  through 
20)  but  has  several  important  differences  and  advantages.  The  Moment  Ordering  Algorithm  has 
as  its  principal  goal  the  discovery  of  a  single  dominant  relationship  in  the  data,  while  the 
Moment  Compression  Algorithm  attempts  to  factor  the  data  into  separable  pieces  or  clusters. 
Two  important  characteristics  that  all  three  of  these  methods  share  is  that  they  openne 
directly  on  the  non-negative  raw  input  matrix  data  and  that  they  reorganize  and  reorder  the 
matrix  data  by  performing  row  and  column  permutations  in  order  to  reveal  obscure  and 


1.  Multiruiate  Auiiyiis  includes  such  mathematical  techniques  u  Regression  Analysis,  Factor  Analysis,  Principal 
Component  Analysis.  Canonical  Analysis,  Ouster  Analysis,  etc. 

2.  CIuMm  Seeking  techniques  ate  those  data  analysis  methods  whkb  soak  to  identify  gtoupt  of  similar  entities. 

Preceding  page  blank 


3 


potentially  informative  data  patterns.  The  output  of  all  these  algorithms,  then,  is  a  new  data 
matrix  with  its  resulting  new  ordering. 

In  Chapter  II,  the  most  important  features  and  characteristics  of  each  of  the  three 
algorithms  will  be  briefly  described.  Then,  in  Chapter  III,  the  major  results  and  conclusions  of 
this  study  will  be  presented.  Part  II  of  this  paper  contains  a  detailed  description  of  the  theory 
and  development  of  the  three  algorithms  along  with  a  number  of  pertinent  examples  which 
illustrate  the  favorable  characteristics  and  general  applicability  of  these  methods. 


II.  DESCRIPTION  AND  OBJECTIVES  OF  THE  THREE  ALGORITHMS 


In  this  chapter  the  three  data  ordering  algorithms  are  briefly  described  and  their 
objectives  are  compared.  More  detailed  description  of  the  theory  and  development  of  these 
algorithms,  along  with  a  number  of  applications,  will  be  found  in  Part  II. 

A.  THE  BOND  ENERGY  ALGORITHM1 


The  Bond  Energy  Algorithm2  is  capable  of  identifying  and  displaying  natural  groups 
and  clusters  that  occur  in  complex  data  matrices.  Moreover,  the  algorithm  is  able  to  uncover 
and  display  the  associations  and  interrelationships  of  these  groups  with  one  another.  These 
tasks  are  accomplished  through  the  use  of  a  numerical  measure  of  how  clustered  or  clumpy3  a 
matrix  is.  The  proposed  measure  of  effectiveness  (ME)  attains  its  maximum  value  when  the 
matrix  assumes  a  very  clumpy  or  aggregated  form.  It  has  beer  found  that  the  structures  and 
relationships  existing  in  data  matrices  more  clearly  exhibit  themselves  when  the  matrices  are 
presented  in  more  aggregated  forms  corresponding  to  larger  MEs. 


The  ME  is  defined  as  follows.  Assume  that  the  matrix  of  relationships  (or  transactions, 
flow,  etc.)  has  dimension  M  by  N  with  non-negative  elements  ay.  The  quantity  ay  is  defined 


aij  J  [ai+  l.j+ai-  l,j+ai.j+  1  +ai.j- 1]' 


From  Fig.  1  it  -an  be  seen  that  ay  is  just  one  half  the  sum  of  the  horizontal  and  vertical 
nearest  neighbors  of  ay.  The  unnormalized  ME  can  now  be  defined  as 


ME  =  L  ay  ay  . 
all  ij 

The  ME  clearly  is  equal  to  the  sum  of  all  the  vertical  and  horizontal  bond  strengths  in  the 
matrix  where  the  strength  of  a  bond  between  two  horizontally  or  vertically  adjacent  elements 

t.  The  theoey  and  development  of  thii  algorithm  are  due  to  t>  W.T,  McCormick,  Jr. 

2.  TnU  algorithm  U  to  called  because  ill  meaiure  of  effective  mm  involve!  product!  of  neareti  neighbor  maim 
element!  that  may  be  likened  to  bond  strength! 

)  A  clumpy  matrix  It  one  who*e  Luge  element!  lie  neat  olhei  boa*  tie  menu,  forming  aggregate*  called  clump*. 

4.  With  the  convention  io^ajo-a^j.a^^j-O. 
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is  defined  as  the  product  of  the  elements.  A  slightly  more  general  form  of  this  ME  is  presented 
in  the  more  detailed  description  in  Part  II  of  this  paper. 

To  obtain  maximum  “dumpiness”  of  the  matrix  it  is  necessary  to  maximize  the  ME 
over  all  row  peosmtations  and  column  permutations  oi  the  matrix,  i.e., 

max 

all  row  perm 
a.  co1  perm 

This  problem  can  be  formulated  equivalently  as  two  quadratic  assignment  problems*  and  its 
optimal  solution  can  be  determined.  However,  this  rigorous  solution  is  quite  time  consuming 
so  a  suboptimai  algorithm  has  been  developed.  The  suboptimal  algorithm  is  a  sequential 
selection  procedure  that  has  proven  to  be  efficient  and  rapid.  The  description  and  details  of 
this  technique  are  contained  u»  Chapter  I  of  Pari  II. 

A  simple  example  will  illustrate  the  sensitivity  of  the  ME  and  the  utility  of  a 
rearrangement  of  the  matrix  data.  Suppose  we  have  a  symmetric  matrix  showing  certain 
associations  or  relationships  between  entities  A,  B,  C  and  D.  The  initial  relationship  matrix  is 
shown  in  Fig.  2a,  where  the  ones  in  the  ij^  elements  of  the  matrix  represent  the  existence  of 
-eiationships  between  entities  i  and  j  and  the  zeros  indicate  the  absence  of  relationships.  It  is 
clear  from  the  definition  of  the  ME  and  the  observation  that  there  are  no  bonds,  that  the  ME  = 
0  for  the  matrix  in  Fig.  2a.  Figures  2b,  2c,  2d,  and  2e  show  progressively  greater  levels  of 
dumpiness  and  their  MEs  are  2,  4,  6,  and  8,  respectively.  Application  of  the  Bond  Energy 
Algorithm  produces  the  ordering  shown  in  Fig.  2e,  where  it  is  dear  that  two  dusters  have  been 

S.  See  Appeadix  K. 
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B.  THE  MOMENT  ORDERING  ALGORITHM6 


The  purpose  of  the  Moment  Ordering  Algorithm  is  to  identify  the  single  dominant 
relationship  in  an  array  of  data,  and  to  reorder  the  rows  and  columns  of  the  array  to  produce  a 
ranking  under  this  dominant  relationship.  That  is  to  say,  the  algorithm  finds  the  principal  axis7 
for  the  data,  and  arranges  both  the  rows  and  the  columns  according  to  the  implicit  underlying 
variable  corresponding  to  the  axis.  The  concept  may  be  made  clearer  by  considering  the  two 
examples  discussed  in  Part  II.  One  example  involves  the  distribution  of  pottery  types  in  a 
group  of  archeological  sites.  The  underlying  variable  is  the  age  of  the  site,  and  the  algorithm 
therefore  produces  a  chronological  ordering  of  the  sites.  The  second  example  involves  the 
voting  patterns  of  a  group  of  Senators.  The  algorithm  determines  that  the  underlying  variable 
is  the  degree  of  liberalism/conservatism,  and  therefore  orders  the  Senators  (and  the  bills  voted 
upon  by  them)  along  a  liberal/conservative  spectrum. 

The  underlying  idea  behind  the  algorithm  is  the  fact  that  if  two  rows  are  very  similar  to 
each  other  their  mean  row  moments  should  be  close  to  each  other  in  value.  The  mean  row 
moment  of  the  ith  row  is  defined  as 

N  /  N 

xi  =  £  i*Hf  /  £ 

j=  1  /  j  =  1 

where  ajj  is  the  ijth  element  in  the  array.  Similarly,  if  two  columns  are  closely  related,  their 
mean  column  moments,  defined  analogously,  should  be  close  in  value.  The  algorithm,  then,  is 
an  attempt  to  use  these  moments  to  rearrange  the  array  so  that  rows  (or  columns)  are  as  near 
as  possible  to  other  similar  rows  (or  columns). 

The  algorithm  begins  by  computing  the  row  moments  for  the  array  in  its  initial  state, 
and  placing  the  rows  in  ascending  order  of  their  moments.  The  column  moments  are  then 
calculated,  and  the  columns  reordered  according  to  their  moments.  This  reordering,  however, 
changes  the  values  of  the  row  moments.  The  row  moments  are  therefore  recalculated  and  the 
rows  reordered.  The  procedure  is  continued,  alternating  between  row  and  column  reorderings, 
until  an  ordering  is  reached  in  which  both  the  rows  and  columns  are  arranged  in  order  of  their 
moments.  Such  a  stable  state  is  considered  a  solution.  The  principal  output  of  the  algorithm  is 
then  the  one-dimensional  ordering  of  the  entities  on  the  array  axes  on  the  basis  of  whatever 
dominant  relationship  may  exist  in  the  data. 


6.  The  Initial  idea  for  thii  algorithm  and  for  thi»  research  paper  is  due  to  Dr.  John  J.  Martin.  The  algorithm  war 
improved  and  developed  by  Dr.  Stephen  B.  Deuttch. 

7.  A  principal  axil  may  be  thought  of  aa  an  “underlying  variable"  by  mean!  of  which  the  explicit  variable*  can  be 
Uited  in  a  ooedimetuionai  ordering. 
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As  an  example  of  how  the  Moment  Ordering  Algorithm  operates  on  a  sample  data 
array,  consider  the  relationship  matrices  given  in  Fig.  3.  When  the  algorithm  is  applied  to  the 
data  array  of  Fig.  3a,  the  new  array  shown  in  Fig.  3b  is  obtained.  Similar  rows  are  now 
adjacent  to  each  other,  and  the  overall  ordering  of  the  rows  reflects  their  placement  along  the 
principal  axis  of  the  array.  Note  the  concentration  of  the  non-zero  elements  along  the  main 
diagonal  of  the  new  array.  This  concentration  is  a  property  of  solutions  found  by  the 
algorithm.  The  details  of  this  method  and  some  examples  which  have  been  successfully 
handled  are  presented  in  Part  II. 
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FIGURE  3.  An  Example  of  the  Moment  Ordering  Algorithm 


C.  THE  MOMENT  COMPRESSION  ALGORITHM8 

The  Moment  Compression  Algorithm  is  designed  to  identify  natural  groups  and  clusters 
of  entities  by  factoring  the  data  reiationshii  mat«x  into  a  number  of  pieces.  The  algorithm 
accomplishes  this  by  finding  the  data  ordering  which  minimizes  a  specific  ME.  The  ME  used  by 
the  Moment  Compression  Algorithm  is  just  the  sum  of  all  the  row  and  column  second 
moments  about  their  respective  means,  that  is 

M  N 

ME»  £  rj  +  £  Cj, 
i=!  j  =  l  1 

where  rj  and  Cj  are  the  i^  row  moment  and  j1^  column  moment.  The  minimization  of  this  ME 
over  all  row  and  column  permutations  has  the  effect  of  compressing  the  data  in  such  a  way  as 
to  force  the  non-zero  matrix  dements  toward  a  block-factored  form. 


1  The  theory  end  development  of  Util  algorithm  ue  due  to  Dr.  Pud  J.  Schweitzer. 
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This  ME  was  devised  because  of  the  observation  that  the  rows  and  columns  of  a  matrix 
in  perfect  block-factored  form,  when  contrasted  with  the  same  matrix  after  row  or  column 
permutations,  have  the  smallest  sum  of  the  moments  of  inertia  about  their  means.  That  is,  any 
row  or  column  permutation  of  a  matrix  in  perfect  block  form  will  “expand”  a  block  and  make 
it  less  dense,  thereby  increasing  the  matrix’s  total  moment  of  inertia.  A  matrix  in  perfect 
block-factored  form  is  shown  in  Fig.  4. 


I M  2-69-6 


FIGURE  4.  Matrix  with  Perfect  Block  Form 


The  problem  of  ME  maximization  can  be  posed  as  two  quadratic  assignment  problems; 
however,  in  practice,  it  has  been  solved  sub-optimally  by  an  iterative  gradient  procedure 
involving  linear  assignment  problems. 

When  the  Moment  Compression  Algorithm  is  applied  to  any  of  the  matrix  orderings  of 
Fig.  2  the  resulting  ordering  is  the  completely  block-factored  form  shown  in  Fig.  2e.  In  this 
special  case  when  the  matrix  is  completely  block  factorable,  the  Bond  Energy  and  the  Moment 
Compression  Algorithms  will  both  produce  block-factored  form. 

D.  CONTRASTS  AMONG  THE  THREE  ALGORITHMS 

In  order  to  understand  better  exactly  how  the  three  algorithms  differ,  it  is  useful  to 
compare  their  objectives  and  their  computational  methods. 

The  single  objective  of  the  Moment  Compression  Algorithm  is  to  identify  groups  or 
clusters  by  rearrangement  of  the  matrix  data.  In  addition  to  sharing  this  objective  the  Bond 
Energy  Algorithm  has  the  additional  objective  of  determining  whether  and  in  what  manner 
these  groups  are  related  to  one  another.9  Computationally,  the  MEs  of  the  two  algorithms 


9.  See  dlicuuioa  before  FI*.  6. 
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differ  substantially  in  that  the  Bond  Energy  ME  depends  on  nearest-neighbor  interactions 
while  the  Moment  Compression  ME  is  completely  global.  A  consequence  of  this  difference  is 
that  the  Bond  Energy  ME  more  adequately  describes  the  topological  properties  of  dumpiness, 
denseness  and  connectedness.  Another  consequence  is  the  greater  computational  ease  in 
optimizing  the  Bond  Energy  ME  by  use  of  a  rapid  sequential  selection  algorithm  which 
exploits  its  nearest  neighbor  dependency. 

The  Moment  Ordering  Algorithm  differs  markedly  from  both  of  the  previous  data 
ordering  methods.  Instead  of  attempting  to  identity  groups,  clusters  or  group  interrelation¬ 
ships,  the  main  objective  of  the  Moment  Ordering  Algorithm  is  to  produce  a  one-dimensional 
ordering  of  entities  along  the  axes  of  the  matrix.  It  accomplished  this  by  finding  a  dominant 
variable  or  principal  axis  along  which  these  entities  can  be  ordered.  Computationally,  like  the 
Moment  Compression  Algorithm,  the  Moment  Ordering  Algorithm  employs  moments  which 
are  global  matrix  measures,  and  thus  it  is  not  as  sensitive  to  local  details  as  the  Bond  Energy 
Algorithm.  Its  principal  computational  difference,  though,  from  the  Bond  Energy  and  the 
Moment  Compression  Algorithms  is  that  it  is  a  completely  heuristic  iterative  technique  that 
docs  not  attempt  to  optimize  any  measure  of  effectiveness. 
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III.  CONCLUSIONS 


The  following  statements  are  the  general  assessments  and  conclusions  regarding  the 
applicability,  overall  usefulness,  and  efficiency  of  the  three  algorithms  developed  for  direct 
analysis  of  multivariate  systems  by  matrix  reordering. 

•  The  Bond  Energy  Algorithm  proved  to  be  the  most  generally  useful  and 
versatile  of  the  three  algorithms  for  treating  certain  problems  of  multivariate 
analysis.  It  is  capable  not  only  of  classifying  and  clustering  data  but  also  of 
successfully  identifying  areas  of  interrelationships  that  exist  among  these 
clusters.  It  has  been  found  to  be  an  efficient  and  general  approach  to  problems 
involving  clusters  and  group  structures. 

•  The  Moment  Ordering  Algorithm  is  an  efficient  technique  for  uncovering  and 
displaying  a  univariate  relationship  inherent  in  the  data.  That  is,  it  is  a  fast  and 
direct  method  for  uncovering  the  principal  axis  of  a  data  structure.  The 
efficiency  of  the  algorithm  was  found  to  be  in  direct  proportion  to  its  ultimate 
success  in  identifying  a  principal  axis.  The  primary  utility  of  this  algorithm  is  in 
determining  a  good  one-dimensional  ordering  of  the  data  rather  than  in 
uncovering  clusters  or  group  interrelationships  in  the  data. 

•  The  Moment  Compression  Algorithm  is  successful  at  identifying  clusters  and 
groups  inherent  in  the  data.  Both  it  and  the  Bond  Energy  Algorithm  will  put  a 
matrix  into  block  form,  if  this  is  possible.  However,  the  Moment  Compression 
Algorithm  is  slower  and  therefore  less  useful  for  large  problems.  Unlike  the 
Bond  Energy  Algorithm,  the  Moment  Compression  Algorithm  cannot  handle 
the  case  of  block-checkerboard1  matrices  arising  from  multilateral  group 
relationships.  Consequently  the  Moment  Compression  Algorithm  is  considered 
inferior  to  the  Bond  Energy  Algorithm  both  with  regard  to  computational 
speed  and  versatility  of  its  measure  of  effectiveness. 


1.  See  Fig.  6. 


Preceding  page  blank 
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I.  THE  BOND  ENERGY  ALGORITHM 

A.  MOTIVATION 

The  motivation  for  the  development  of  the  Bond  Energy  Algorithm  was  to  be  able  to 
treat  a  broader  class  of  problems  than  that  normally  found  in  cluster  analysis  applications.  In 
addition,  it  was  desired  to  operate  directly  on  and  manipulate  the  original  data  without 
creating  or  losing  information.  The  object  is  not  only  to  classify  and  group  similar  entities  but 
also  to  determine  how  and  by  what  means  these  groups  are  interrelated.  This  can  be  illustrated 
by  considering  a  symmetric  binary  (0-1)  relationship  matrix  between  N  entities.  If  the  N 
entities  can  be  separated  into,  say,  four  unique  groups  (unique  meaning  that  the  entities  in  one 
group  are  related  only  among  themselves  and  not  with  any  entities  outside  their  own  group), 
then  many  of  the  techniques  of  cluster  analysis  are  applicable.  In  this  case  it  is  possible  to 
reorder  the  rows  and  columns  of  the  input  data  matrix  to  obtain  the  form  given  in  Fig.  5. 

-  CN  Z 


n«i2^; 


FIGURE  5.  Relationship  Matrix  Showing  4  Unique  Groups 

However,  if  the  entities  are  not  completely  factorable  into  unique  groups  then  it  is  often 
desirable  to  identify  not  only  the  principal  groups  but  also  their  significant  areas  of  relation¬ 
ship.  In  other  words,  it  might  be  desirable  to  rearrange  the  data  matrix  to  obtain  a  checker¬ 
board  pattern  if  it  is  possible.  This  type  of  pattern  is  shown  in  Fig.  6,  where  the  off-diagonal 
blocks  of  large  Xs  represent  data  clumps  containing  a  sizable  percentage  of  non-zero  entries, 
thus  indicating  partial  or  total  intergroup  relationships. 

Preceding  pap  blank 
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FIGURE  6.  Relationship  Matrix  with  Block-Checkerboard  Form 

The  essential  question  is,  given  a  matrix  whore  the  data  are  presented  in  an  arbitrary  manner, 
how  can  the  rows  and  columns  of  a  matrix  be  simply  rearranged  to  obtain  as  “clumpy”  a 
matrix  form  as  possible. 

B.  THE  MEASURE  OF  EFFECTIVENESS 

1 .  Definition  and  Interpretations 

In  order  to  analytically  determine  the  “dumpiness”  of  a  particular  matrix,  it  was 
necessary  to  develop  some  measure  of  effectiveness  (ME) 1  of  how  any  subsequently  proposed 
algorithm  would  progress.  This  ME  must  be  sensitive  to  and  depend  on  local  dumpiness  while 
also  characterizing  the  dumpiness  of  the  entire  matrix.  The  essential  idea  behind  the  ME, 
which  fulfills  this  requirement,  came  from  likening  the  situation  to  that  of  the  saturation  of 
bonds  in  the  nucleus  of  an  atom.  That  is,  when  the  nucleons  are  clumped  together  there  is 
total  bond  saturation  in  the  interior  of  the  nucleus  while  the  bonds  of  the  nucleons  near  the 
surface  are  unsaturated.  The  intent  was  to  find  an  ME  which  when  maximized,  resulted  in  as 
few  unattached  or  unbonded  matrix  elements  as  possible.  The  bond  strength  between  two 
adjacent  matrix  elements  is  defined  as  the  l/k^  power  of  the  product  of  the  matrix  elements. 
Maximization  of  the  ME  will  maximize  the  sum  of  all  the  bond  strengths,  and  therefore  clump 
together  the  larger  non-zero  matrix  elements.  Another  physical  phenomenon  that  may  be 
likened  to  this  situation  is  that  of  water  beads  on  a  glass.  The  beads  tend  to  aggregate  into 


1.  A  mote  complete  dtomlon  of  StE»  aty  b«  found  la  Appendix  G. 


larger  clumps  in  order  to  minimize  the  surface  energy.  The  ME  can  be  defined,  then,  as  just  the 
sum  of  all  the  bond  strengths  in  a  matrix.  Thus 


MEk=  E  a.1/1'  ,  a.. 
k  all  ij  »  k  » 


where2 


kaij 


=  ±  fa1/k 

2  aij+l 


.  1/k  ,  1/k  . 

i. . . ,  +  a. .  ,  +  a. . ,  .  +  a 
U+l  y-1  i+l  j 


1/k' 

i-lj_ 


and  k  is  a  weighting  constant,  which  is  usually  set  equal  to  2.  The  ME  may  be  interpreted 
mathematically  as  the  sum  of  the  scalar  products  (or  projections  on  one  another)  of  all  the 
contiguous  row  vectors3  plus  the  sum  of  the  scalar  products  of  all  of  the  contiguous  column 
vectors.3 


2.  Normalization  of  the  ME 

The  ME  defined  above  can  be  normalized  so  that  its  value  varies  between  0  and  1 .  This 
normalized  measure  of  effectiveness  (NME)  is  defined  as 


NMEk  =  i  E  a!/ka.. 
K  ^  ij  k  y 


all  ij 


and 


0  <  NMEk  <  1 


where 

S  is  the  normalization  constant  defined  as 
S  -  2  E  af . 

all  ij  u 

S  can  be  interpreted  mathematically  as  the  sum  of  the  squares  of  the  L->  norms4  of  all  the  row 
and  column  vectors.  The  advantage  in  having  a  normalized  ME  is  that  it  is  easier  to  determine 
how  much  improvement  in  the  dumpiness  of  a  matrix  has  been  achieved  since  it  is  a  measure 

2.  Apln.^"  *y,  »  »M*1J  *  \N*t 

J.  The  I**1  row  vector  U  oomptUcd  of  the  etemeati  , ...  in  the  i1^  low  of  the  matrix.  A  column  vector 

l*  defined  anaiofourly. 

4.  The  Lj  vector  norm  ii  defined  n  / - 

H  :/e 

V  l»  I 

where  M  U  the  dimetuion  of  the  vector  tpaee.  The  fact  that  the  NME^  U  property  normatited  foUowi  from  a  hade  inequality 
for  a  tvormed  tpece  |  ♦  |bp3*2(a.^). 
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of  the  amount  of  bond  saturation.  For  instance,  if  the  NME  of  the  reordered  data  matrix 
equals  0.6,  whereas  the  NME  of  the  initial  data  matrix  equals  0.2,  then  it  can  be  concluded 
that  there  does  exist  a  good  deal  of  inherent  group  structure  and  interrelationship  that  was  not 
initially  evident.  Moreover,  the  final  N’iE  gives  an  absolute  measure  of  the  existence  of  the 
clusters  that  we  have  sought  to  uncover. 

3.  Advantages  of  the  ME 

The  ME  proposed  above  has  some  very  important  theoretical  and  computational 
advantages  which  will  be  enumerated  here. 

•The  NME  can  be  used  for  matrices  of  any  size  or  shape.  In  addition,  symmetry 
of  the  matrix  is  not  required.  The  only  restriction  is  that  the  matrix  elements 
be  non-negative,  real  numbers. 

•  Since  the  vertical  (horizontal)  bonds  are  unaffected  by  interchange  of  the 
columns  (rows),  the  ME  decomposes  into  two  parts;  one  (sum  of  the  vertical 
bonds)  dependent  only  on  row  permutations,  and  the  other  (sum  of  the 
horizontal  bonds)  dependent  only  on  column  permutations.  Consequently 
optimization  of  the  ME  can  be  achieved  in  exactly  two  passes,  one  finding  the 
optimal  column  permutation,  the  other  finding  the  optimal  row  permutation. 

•  These  two  passes  can  be  carried  out  completely  independently  of  each  other. 
In  particular,  it  is  not  necessary  to  alternate  between  row  and  column  permu¬ 
tations,  as  in  the  Moment  Ordering  Algorithm,  thus  eliminating  the  possibility 
of  any  cycling5  of  the  solution. 

•  Since  the  contribution  to  the  ME  from  any  column  (or  row)  is  only  affected  by 
the  two  adjacent  columns  (or  rows),  the  optimization  lends  itself  very  well  to  a 
multistage  sequential  selection  process. 

•lire  Bond  Energy  ME  optimization  does  not  require  any  prior  prejudices,  such 
as  forcing  the  data  into  clumps  along  the  diagonal  or  forcing  the  data  into 
block-diagonal  form.  The  representation  of  the  data  that  is  sought  is  a  tight 
clumped  form  and  so  the  maximization  of  the  ME  might  very  well  allow  the 
possibility  of  far  outlying  elements  in  order  to  achieve  globally  higher  degree  of 
compactness.  Tlus  feature  is  partieuk  rly  important  in  the  case  of  multilateral 
relations  between  groups  of  entities  where  i!  is  clearly  not  possible  to  obtain  a 
block-diagonal  form. 

3.  TWi  phenomena »  oocuii  whan  (be  totulkw  |eu  aUcratiely  belie*  then  wo»«e. 
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•The  1/k  power6  of  ay  appearing  in  the  expression  for  the  ME  allows  any  desired 
weighting  of  the  larger  matrix  elements. 

C.  THE  SOLUTION:  MAXIMIZATION  OF  THE  ME 

1.  The  Fxact  Solution 

The  problem  to  be  solved  as  implied  earlier  is  to  maximize  the  ME  over  all  row  and 
column  permutations.  That  is, 

*“  ji  E  a'/k  fa'*  +  al/k 

***  *2  >.j  Jr(i),0(j)  [*(0,00  +  1)  *(0,  00  -  1) 

+  al/k  + 

»(i  +  1),  00)  ff(i-l),0(j) 

where  *  =  {jt(1),  ir  (2), . .  .  ,  ir(M)}  and  0  =  {  0(1),  0(2),  . .  .  ,  0(N) } 

are  the  row  and  column  permutations.  This  can  be  thought  of  physically  as  maximizing  the 
sum  of  all  the  bond  energies  and  mathematically  as  maximizing  the  sum  of  all  the  scalar 
*■  products  of  contiguous  row  vectors  and  column  vectors.  This  maximization  problem  can  be 
stated  equivalently  as  two  quadratic  assignment  problems  (the  reader  is  referred  to  Appendix 
A  for  the  detailed  formalism).  The  first  seeks  a  permutation  of  the  columns  of  [  ay )  which 
maximizes  the  row  bond  energy,  the  other  seeks  a  permutation  of  the  rows  of  [ ay]  which 
maximizes  the  column  bond  energy.  These  optimizations  may  be  viewed  as  two  clustering 
procedures,  one  which  reorders  the  rows  on  the  basis  of  their  similarity  (similarity  being 
measured  by  the  scalar  product  of  the  two  rows),  the  other  reordering  the  columns.  Reas¬ 
sembling  the  matrix  after  both  clusterings  produces  the  dense  blocks  shown  in  Fig.  6. 
Although  quadratic  assignment  problems  can  be  solved  exactly  as  well  as  approximately  (for 
exact  and  approximate  methods  see  the  references  listed  in  Appendix  A),  the  solution  of  this 
problem  requires  a  large  amount  of  computer  time  in  either  case.  Our  own  approximate 
sequential  selection  algorithm  has  been  developed  which  takes  advantage  of  the  nearest 
neighbor  properties  of  the  measure  of  effectiveness. 

2.  Approximate  Solution 

a.  Description  of  the  Sequential  Selection  Algorithm  The  suboptimal  algoritlun  which 
has  been  actually  used  to  determine  local  optima  of  the  ME  is  as  follows: 

~  d.  rhe  twulttvUy  of  the  ME  to  the  vtiue  of  -  U4taeu*»»iin  Appendix  G. 
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(1)  Compute  and  store  the  scalar  products  of  each  row  with  every  other 
row  and  each  column  with  every  other  column. 

(2)  Select  any  column  to  begin  the  selection  process.  Set  i=l . 

(3)  Next,  try  each  of  the  remaining  N-l  columns  placed  alongside  the  first 
column  and  compare  its  contribution1 2 3 4 5 6 7  to  the  horizontal  bond  ME. 

(4)  Place  alongside  the  first  column  that  particular  column  which  gives  the 
largest  contribution  to  the  ME. 

(5)  Continue  the  process  at  the  i^  step  by  comparing  the  contribution  to 
the  ME  by  placing  each  of  the  N-i  remaining  columns  in  each  of  the  i+1 
possible  positions,®  and  putting  the  one  which  gives  the  largest  con¬ 
tribution  to  the  ME  in  its  proper  place. 

(6)  After  the  process  is  completed  by  placing  the  last  remaining  column  in 
its  “best”  place,  then  the  entire  procedure  (items  2  through  5)  is 
repeated  on  the  rows.  It  is,  however,  not  necessary  to  repeat  the 
procedure  on  the  rows  if  the  initial  input  matrix  is  symmetric  since  the 
final  resulting  row  order  will  be  identical  with  the  column  ordering, 
yielding  a  symmetric  matrix. 

b.  Advantages  of  the  Algorithm.  The  algorithm  described  above  has  several  attractive 
advantages  which  are  noted  here. 

(1)  Since  the  algorithm  is  finite  and  non-iterative,  there  are  no  convergence 
problems. 

(2)  The  algorithm  will  always  reduce  a  matrix  to  pure  block  form  if  it  is 
possible  to  obtain  this  form  by  row  and  column  permutations  (see 
Appendix  B  for  proof). 

(3)  The  solution  obtained  from  the  algorithm  is  independent  of  the  input 
order  of  the  rows  (or  columns)  but  is  only  dependent  on  the  initial  row 
(or  column)  chosen  to  start  the  sequential  selection  process. 

(4)  The  results  of  the  algorithm  are  very  insensitive  to  the  starting  point 
(i.e.,  starting  row  or  column),  hence  any  solution  is  a  “good”  one  (see 
Tabic  2). 

(5)  The  computation  time  for  the  algorithm  depends  only  on  the  size  of 
the  matrix  ano  not  on  its  elements. 

(6)  The  algorithm  uses  no  thresholds  or  filtering  during  its  operation  which 
can  alter  its  course  and  affect  the  final  result. 

(7)  Only  the  raw  input  data  matrix  is  used  to  determine  the  new  row  and 
column  orderings. 

7.  The  coauiOutloa  U  )u»i  Uu  4oi  product  of  the  chown  column  veste?  with  the  flut  column  veetot. 

8.  The  HI  podUotu  we  to  the  left  end  rtjhs  of  the  1  eotumm  already  pitced. 


3.  An  Example 


A  simple  example  taken  from  Principles  of  Numerical  Taxonomy  (Ref.  1 5)  will 
illustrate  how  the  algorithm  can  identify  the  clusters  and  their  interrelationships.  The 
similarity  matrix  of  Fig.  7  is  given  where  a  numerical  value  of  5  in  element  i  j  indicates  a  high 
degree  of  similarity  between  entity  i  and  entity  j,  and  0  indicates  no  similarity. 

A  B  C  D  E  F  G  H  I  J 

A  5  4  1  0  4  1  1  0  3  1 

B  4  5  0  1  3  1  1  0  4  1 

C1050133012 
D0105000401 
E  4  3  1  0  5  1  0  0  4  1 

F  l  1  3  0  1  5  3  0  1  3 

G  1  130035012 

H0004000510 
I  3  4  1  0  4  1  1  1  5  1 

J  1  1  2  1  1  3  2  0  1  5 

I M  2-69-9 

FIGURE  7.  Initial  Non-Binary  Similarity  Matrix 
By  applying  the  algorithm  described  above,  a  new  axis  ordering  and  a  new  matrix  are  obtained 
and  are  shown  in  Fig.  8. 
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FIGURE  8.  Reordered  Non-Binary  Similarity  Matrix 


It  is  easy  to  identify  three  major  diagonal  blocks  of  large  numbers  representing  three  dusters 
or  groups  of  entities.  H  and  O  constitute  the  first  group,  B,  A,  E  and  1  the  second  group,  and  J, 
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F,  G  and  C  the  third  group.  From  the  grouping  of  the  smaller  off-diagonal  elements  it  is 
evident  that  there  is  some  weak  relationship  between  the  second  and  third  groups  but 
essentially  no  relationship  between  the  first  group  and  either  of  the  other  two.  It  is  also  quite 
apparent  from  this  example  that  this  new  form  for  the  matrix  data  conveys  more  information 
concerning  the  group  structure  and  relationships  than  does  the  original  matrix  form. 

D.  OPERATION  OF  THE  ALGORITHM 

1.  Computing  Time  Requirements 

If  the  original  data  matrix  is  of  dimension  M  by  N,  then  the  total  number  of  arithmetic 
operations  necessary  to  perform  all  the  initial  row  and  column  dot  products  is  just: 

M- 1  N- 1 

Operations  =  N  E  i  +  M  V  j  or, 

i-  1  j  =  1 

Operations  =  N  *  -  I)  +  ^  .  N][N  -  1)  or 

2  2 

for  large  M  and  N, 

Operations  =  M^N  +  N2M 

At  step  i  of  the  algorithm,  it  is  necessary  to  compare  the  contribution  of  the  ME  of  all  the 
remaining  N-i  unplaced  columns  in  the  i+1  possible  positions,  thus  the  total  number  of  column 
comparisons  equals 

N  -  1  N  - 1 

E  (i+  1 )  (N  -  i)  =  E  iN-i2  +  N-i 

i  = 1  i=l 

N3 

for  large  N. 

6 

Similarly  it  requires  approximately  M3/6  comparisons  for  the  rows.  Thus  for  a  square  matrix, 
the  computation  time  of  the  algorithm  goes  as  N3.  This  theoretical  variation  in  the  computing 
time  has  been  borne  out  experimentally  as  can  be  seen  in  Table  1.  The  computing  time9  in 
seconds  is  given  for  various  size  problems  (times  given  are  for  a  single  startup  point). 


9.  Ob  CDC  1604  eompuUt  uiisf  Uoivasily  of  KIobcjoU  compiler. 
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Table  1.  BOND  ENERGY  ALGORITHM  COMPUTATION  TIME 


N 

M 

Time 

21 

21 

1 1  sec 

29 

29 

23  sec 

48 

48 

1 24  sec 

From  this  data  a  scaling  law  can  be  derived  which  gives  the  required  computation  time 
in  seconds  for  a  given  size  square  matrix,  for  a  single  starting  point. 

Computation  Time  (sec)  =  0.001 2 

2.  Ties 

It  occasionally  happens  that  ties  occur  during  the  course  of  the  sequential  selection 
algorithm.  Ties  between  rows  and  columns  can  occur  in  the  following  ways: 

(1)  Tic  arising  from  putting  the  same  as  yet  unplaced  column  (or  row)  in 
two  or  more  possible  positions. 

(2)  Tic  arising  from  putting  different  as  yet  unplaced  columns  (or  rows)  in 
the  same  positions. 

(3)  Ties  arising  from  putting  different  as  yet  unplaced  columns  (or  rows)  in 
different  possible  positions. 

We  have  no  present  criterion  for  deciding  how  to  break  ties  arising  from  condition  ( 1 ),  nor  is  it 
known  whether  there  is  reason  to  select  one  alternative  over  the  others.  Ties  arising  from 
conditions  (?.)  and  (3)  are  broken  by  selecting  the  unplaced  row  or  column  which  has  the 
shortest  length  in  the  L2  norm.10  Thinking  in  terms  of  the  ME  mathematically,  if  we  can 
obtain  the  same  scalar  products  or  projections  witn  two  vectors,  then  the  shorter  should  be 
used  rather  than  the  larger  one.  This  tie-breaking  mechanism  has  been  found  to  work 
satisfactorily  in  that  it  leads  to  informative  final  data  arrangements. 

3.  Effect  of  Starting  Point 

Although  the  results  of  the  algorithm  do  not  depend  on  the  order  in  which  the  rows 
and  columns  are  considered,  there  is  a  difference  in  the  final  results  depending  on  which  row 
or  column  is  selected  to  start  the  multistage  decision  process.  In  the  example  presented  in  Figs. 
7  and  8,  the  problem  was  started  10  times,  beginning  once  with  each  column.  Table  2  gives  the 
frequency  of  occurrence  and  final  ME  for  each  distinct  solution. 

t0.  TttU  U  juU  Uw  square  root  of  the  sum  of  the  tquuiet  of  eti  the  eteuMniv  of  the  vector. 
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Table  2.  FREQUENCY  DISTRIB  UTION  OF  ME 


Solution  No. 

Frequency 

ME 

1 

3 

419 

2 

3 

419 

3 

1 

414 

4 

2 

414 

5 

1 

412 

Several  significant  facts  may  be  noted  from  these  results.  First  of  all,  the  solutions  with  the 
highest  ME,  419,  which  are  believed  to  be  the  globally  optimum  solutions,  occur  60  percent  of 
the  time.  The  difference  between  the  best  and  worst  solution  is  only  7  out  of  over  400,  or  less 
than  2  percent.  A  noteworthy  point  here  is  that  the  final  “solution”  (ME)  depends  very 
weakly  on  the  starting  point  and  even  the  worst  “solution”  is  not  very  far  from  the  optimal 
solution.  With  regard  to  the  final  group  structure,  it  has  been  found  that  the  various  near 
optimal  so'utions  do  not  differ  significantly  from  the  optimal  solution.  The  various  solutions 
are  due  to  the  rearrangement  of  the  entities  within  a  cluster  group  and  the  reordering  of  the 
groups  themselves.  These  results  have  been  confirmed  by  experimentation  on  significantly 
larger  matrices. 

4.  Formatting  Data 

The  input  format  for  the  data  can  be  in  any  matrix  form.  This  means  that  the  Bond 
Energy  Algorithm  permits  analysis  of  the  raw  data  without  forming  a  similarity  matrix.*  *  For 
example,  suppose  we  have  an  object-attribute  matrix  and  we  desire  to  find  out  which  objects 
are  similar.  The  advantage  of  performing  the  grouping  directly  upon  the  object-attribute 
matrix,  rather  than  upon  the  similarity  matrix,  is  that  it  is  now  possible  to  determine  which 
attributes  characterize  a  particular  group  of  objects  (see  example  4). 

E  APPLICATIONS 

Several  applications  of  this  method  have  already  been  attempted  and  others  have  been 
suggested.  It  appears  that  the  algorithm  is  applicable  to  a  wide  class  of  problems,  a  number  of 
wnich  will  be  enumerated  here. 

( 1 )  Identification  of  natural  groups  and  subgroups  within  data. 

(2)  Identification  of  relationships  and  dependencies  between  groups. 

(3)  Relationships  of  groups  of  attributes  to  groups  of  objects. 

(4)  Examining  influence  relationships  and  structures  via  nonsymraetrical 
data  matrices. 

II.  A  Xmllutty  rntfix  U  a tymaeute  mUix  wfeoie  ijlth  element  U a  meuwe  oT  the  Umikrity  of  entity  i  to  entity  j. 
Appiyti*  the  Bond  En«*y  Afcoeithm  to  a  Omikriiy  iwautt  ideaiiitei  (u  the  diatrsoai  hiociu)  the  mein  of 

entitkt  end  (u  the  ofMiaaooal  dump*)  the  iateqpoup  lelaiiooihipa. 
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(5)  Analysis  of  hierarchical  clustering  and  grouping  via  quantified  numeri¬ 
cal  relationships. 

(6)  Factoring  of  large  linear  assignment  problems  (Ref.  32). 

(7)  Factoring  of  large  management  problems  to  identify  optimal  subtasks. 

(8)  Clustering  of  correlation  matrices. 

(9)  Solution  of  traveling  salesma  n  problems  ( Ref.  3 1 ). 

( 1 0)  Unscrambling  flow  graphs  and  network  relationships. 

F.  EXAMPLES 


A  number  of  examples  are  presented  in  the  following  paragraphs  to  illustrate  the 
operation  and  the  potential  application  of  the  Bond  Energy  Algorithm.  It  should  be  clearly 
understood  that  the  algorithm  operates  on  matrices  that  contain  “hard”  numerical  entries  and 
therefore  considers  each  data  matrix  to  be  an  exact  representation  of  the  relationships 
involved.  We  feel,  nevertheless,  that  the  algorithm  has  application  for  problems  involving 
“soft”  data  (Airport  example)  as  well  as  “hard”  data  (Hindi  consonant  example),  as  long  as 
proper  care  is  taken  to  judiciously  weigh  the  results  subject  to  the  degree  of  validity  of  the 
input  information. 


1.  Example  1 

Bonner  (Ref.  3)  has  presented  several  clustering  techniques  which  uncover  group 
structure  in  matrix  data.  For  this  example,  the  Bond  Energy  Algorithm  is  applied  in  several 
different  ways  to  illustrate  its  advantages  and  directness  for  gathering  similar  data  into  clusters. 
The  objects  which  will  be  clustered  are  defined  by  a  set  of  attributes  which  characterize  them. 


Bonner  presents  a  binary  description  of  an  object  set  as  an  object-attribute  matrix 
which  is  shown  in  Fig.  9.  ATTR,  BUTE  NUMBER 
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FIGURE?.  Initial  Binaiy  Object-Attribute  Matrix 
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He  then  proceeds  to  form  a  similarity  matrix  P,  where  the  Py  are  defined  as 


CH  + 


Sl_ 

cii'  ^ 


and  Cy  is  the  number  of  attributes  which  are  “one”  for  both  object  i  and  object  j.  The 
similarity  matrix  corresponding  to  Fig.  9  is  shown  in  Fig.  10.  A  threshold  T=0.45  is  then  used 
to  convert  the  fractional  similarity  matrix  of  Fig.  10  to  a  binary  similarity  matrix  by  setting 
those  matrix  elements  to  one  whose  values  are  greater  than  0.45  and  the  rest  equal  to  zero. 

i 

This  similarity  matrix  is  shown  in  Fig.  1 1.  Bonner  then  uses  this  similarity  matrix  as  a  starting 
point  for  several  clustering  techniques. 
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FIGURE  10.  Initial  Fractional  Object  Similarity  Matrix 
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FIGURE  II.  Initial  Binary  Similarity  Matrix 


The  Bond  Energy  Algorithm  has  several  advantages  over  Bonner’s  technique.  First,  it  is 
able  to  operate  directly  on  the  object-attribute  matrix  without  forming  a  similarity  matrix  thus 
permitting  it  to  identify  those  particular  attributes  that  characterize  objects  in  the  same 
cluster.  Second,  the  application  of  the  algorithm  will  never  result  in  a  loss  of  information  since 
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the  data  are  only  rearranged.  Finally,  even  using  a  similarity  matrix  the  algorithm  can  produce 
a  reordering  which  not  only  displays  the  clusters  but  also  their  strengths  and  relationships. 
These  advantages  will  be  demonstrated  by  successive  application  of  the  Bond  Energy  Algo¬ 
rithm  to  the  matrices  of  Figs.  9,  10  and  1 1. 


When  the  object-attribute  matrix  of  Fig.  9  is  rearranged  by  the  Bond  Energy 
Algorithm,  the  new  data  matrix  of  Fig.  1 2  is  obtained.  When  rectangles  arc  constructed  around 
solid  blocks  of  Is  in  two  or  more  rows  and  columns,  it  can  be  seen  that  the  objects  fall  into  4 
“core”  clusters:  3,6  and  2,1,5,  and  4,7,  and  8.  it  is  also  observed  that  attributes  3  and  5  arc  the 
essential  characterizing  attributes  of  the  3,6  object  cluster,  attributes  4  and  1  arc  the 
characterizing  ones  for  the  cluster  containing  objects  2,1  and  5,  and  attributes  2  and  6 
characterize  the  cluster  containing  objects  4  and  7. 
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FIGURE  12,  Reordered  Binary  Object-Attribute  Matrix 


When  the  Bond  Energy  Algorithm  is  applied  to  the  fractional  object  similarity  matrix 
of  Fig.  10,  a  new  ordering  is  obtained.  In  this  new  ordering  in  Fig.  13, 


FIGURE  13.  Reordered  Fractionol  Object  Similarity  Matrix 
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only  the  larger  elements  (i.e.,  1/2  or  greater)  are  shown  so  that  the  clusters  can  be  more  easily 
identified.  Again,  it  is  possible  to  identify  the  clusters  and  how  they  interrelate.  Objects  3  and 
6  form  a  very  tight  independent  cluster.  Objects  5,1,2  form  another  tight  cluster,  although 
there  is  a  non-trivial  relationship  between  object  2  and  objects  7  and  8.  Objects  4  and  7  form 
another  cluster  that  is  somewhat  related  to  objects  2  and  8.  Thus,  visually,  this  form  of  data 
presentation  is  helpful  and  its  computational  requirements  are  very  small  (less  than  one 
second). 

When  the  Bond  Energy  Algorithm  is  applied  to  the  binary  similarity  matrix  of  Fig.  1 1, 
the  result  is  given  in  Fig.  14.  It  is  quite  apparent  that  these  results  illustrate  the  same 
relationships  and  clusters  as  those  shown  in  Fig.  13,  but  are  inferior  since  the  strengths  of  the 
relationships  are  not  shown.  This  illustrates  that  while  Bonner’s  filtering  technique  leads  to  the 
uncovering  of  major  clusters,  it  also  loses  information  present  in  the  original  data  matrix. 

OBJECT  NUMBER 


OBJECT  NUMBER 


FIGURE  14.  Reordered  Binary  Similarity  Matrix 
2.  Marketing  Techniques  and  Applications 

In  displaying  the  data  relationsliips  in  this  example,  it  is  found  that  the  application  of 
the  Bond  Energy  Algorithm  reveals  several  latent  group  associations  and  significantly  enhances 
the  quality  of  the  presentation  of  the  data. 

Figure  15  contains  a  matrix  showing  which  Marketing  Techniques  are  used  for 
particular  Marketing  Applications. 11  By  application  of  the  algorithm  it  is  possible  to  reorder 
or  relist  the  marketing  applications  on  the  one  axis  and  the  marketing  techniques  on  the  other 

12.  The  data  tor  this  example  were  taken  from  the  SaptembwOctobai  1969 1 saw  of  the  Harvard  Bualnear  Review 
from  “Technique*  in  Marketing  Rt  March"  by  J.F,  Duh  and  C.  Bert  won. 
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FIGURE  15.  InIHi 
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axis,  while  preserving  all  the  data  relationships  contained  in  the  original  matrix.  Fig.  16 
contains  the  reordered  matrix  produced  by  the  Bond  Energy  Algorithm.  With  the  data  in  this 
new  matrix  form,  it  is  now  possible  to  identify  three  major  clusters  or  clumps  of  data.  It  is 
believed  that  in  this  new  form  it  is  possible  to  uncover  useful  information  that  was  not  obvious 
from  the  original  matrix. 

First,  the  algorithm  groups  marketing  analysis  techniques  that  are  used  for  the  same 
applications  and  also  it  groups  marketing  applications  that  utilize  the  same  marketing  tech¬ 
niques.  This  has  the  effect  of  putting  similar  marketing  techniques  near  one  another  on  the 
vertical  axis  and  putting  similar  applications  together  on  the  horizontal  axis.  It  is  postulated 
that  the  clumps  provide,  for  one  thing,  a  basis  for  efficient  assignment  of  responsibilities  to 
analysts  and  their  supervisors,  and  for  another,  by  exception,  a  basis  for  deciding  upon  the 
relative  merits  of  “techniques"  specialists  and  “applications"  specialists. 

Second,  if  it  is  possible  to  factor  a  matrix  completely  so  that  it  is  apparent  that  there  is 
a  unique  relationship  between  a  certain  group  of  marketing  techniques  and  a  certain  group  of 
marketing  applications,  then  the  algorithm  will  accomplish  this.  In  this  example,  this  has  been 
partially  done  by  identifying  three  more  or  less  independent  clumps  in  Fig.  16.  In  particular,  as 
was  noted  by  the  authors,  PERT  and  CPM  are  similar  in  concept  and  hence  they  occur 
together  in  the  same  clump.  Also,  as  noted  in  the  article,  risk  analysis  is  often  used  in 
conjunction  with  the  method  of  decision  trees.  Here  again,  these  marketing  techniques  are 
contiguous  in  the  new  ordering.  On  the  other  axis  it  is  found  that  similar  “marketing 
applications"  are  grouped  together.  For  example,  Product  planning,  R&D  planning.  Venture 
planning  and  Product-line  analysis  all  in/olve  planning  of  some  sort  and  occur  in  the  same 
clump  because  they  utilize  common  “marketing  techniques”  for  planning,  such  as  PERT,  CPM, 
etc. 


Another  possible  way  in  which  the  clumped  matrix  of  Fig.  16  can  be  useful  is  to 
suggest  possible  unexploited  application  of  techniques  to  marketing  applications  to  which  they 
have  not  already  been  applied.  These  could  be  identified  by  looking  for  conspicuous  holes 
within  the  dumps  or  omissions  on  the  borders  of  the  dumps. 

Thus,  it  appears  in  this  example,  that  with  proper  arrangement  of  binary  (yes-no)  or 
quantified  data  given  in  matrix  representation,  that  the  amount  of  information  conveyed  can 
be  significantly  enhanced  to  such  an  extent  that  it  is  undesirable  to  present  it  in  other  than 
clustered  form. 


Preceding  page  blank 


3.  Coordinating  Airport  Design13 


A  practical  way  to  design  an  airport  is  to  factor  the  problem  into  a  number  of  smaller 
pieces.  If  the  subproblems  can  be  solved  separately  and  then  adjusted  so  as  to  remain  valid  in 
the  context  of  the  original  problem,  then  the  task  is  completed.  It  is,  however,  necessary  to 
determine  the  best  way  to  factor  the  big  problem  into  more  manageable  pieces. 

A  numerical  example  will  illustrate  the  applicability  of  the  Bond  Energy  Algorithm  to 
the  problem.  The  objective  is  to  exploit  the  structure  of  an  airport  problem  in  such  a  way  as  to 
identify  two  things: 

•The  “natural”  subproblems 

•  The  necessary  coordination  between  subproblems. 

The  ultimate  accomplishment  would  be  to  factor  the  problem  into  small,  completely 
independent  subproblems.  But  given  that  complete  independence  is  impossible,  the  next  best 
thing  is  to  minimize  the  intergroup  dependencies  by  identifying  the  optimal  way  to  subdivide 
the  problem. 

The  first  step  is  to  describe  the  airport  problem  in  terms  of  a  set  of  variables  and  their 
interrelations.  A  partial  list  of  exogenous  and  control  variables  is  shown  in  Table  3. 

The  exogenous  variables  describe  those  factors  mostly  dictated  by  the  environment 
while  the  control  variables  apply  to  those  factors  primarily  under  control  of  an  airport  planner. 
Let  Xj  be  the  i1*1  exogenous  variable  and  let  Dj  be  the  i1^  control  variable.  The  Xj’s  may  be 
thought  of  as  input  data  and  Dj’s  as  the  design  decisions.  Given  a  set  of  values  for  the  Xj’s,  it  is 
assumed  that  there  exists  some  way  of  measuring  the  performance  of  an  airport  design  based 
on  some  criteria.  The  details  of  the  performance  function  are  not  needed:  just  a  few  basic 
characteristics.  Let  P  be  the  measure  of  performance  and  let  F  be  the  function  that  measures 
performance.  Geariy,  P  is  a  function  of  the  Dj’s.  hence 

P=F(D1,D2...>D27), 

F  will,  in  general,  depend  on  the  Xj’s;  however,  the  discussion  will  be  limited  to  a  specific  set 
of  values  for  the  Xj’s.  The  design  problem  involves  selecting  values  for  the  Dj’s  that  maximize 

t ).  The  analysis  wul  -he  data  for  this  application  are  due  to  Mr.  T.W.  White  of  the  Institute  for  Defense  Analyses. 


Table  3.  AIRPORT  VARIABLES 


Kxogenous  Variables 

!  Total  air  travel  demand 

2.  Originating  passengers 

3.  Transferring  passengers 

4.  Terminating  passengers 

5.  Greeters  and  well-wishers 

6.  Access  ground  transportation  mode  for  passengers 

7.  Egress  ground  transportation  mode  for  passengers 

8.  Airport  employees 

9.  Taxis  and  cars  that  do  not  park 

10.  Cars  whose  drivers  park  and  fly 

1 1 .  Renta!  cars  going  to  the  airport 

1 2.  Rental  cars  driven  from  the  airport 

13.  Bus  and  limousine 

14.  Employee  access  transportation  mode 

1 5.  Passenger  trip  duration 

16.  Aircraft  turn  around  time  on  apron 

1 7.  Mix  of  aircraft  by  capacity 

IS.  Gate  schedule:  aircraft  arrivals  and  departures 

19.  Origin/destination  pattern  for  baggage  at  airport 

20.  Air  cargo  demand 

21.  Runway  demand 

Control  Variables 

1 .  Passenger  check-in 

2.  Baggage  check-in 

3.  Baggage  claim 

4.  Baggage  moving  system 

5.  Intra-airport  transportation  system 

6.  Cargo  terminal 

7.  Close-in  parking  lots 

8.  Remote  parking  lots 

9.  Main  access  roads  to  and  from  airport 

10.  Circulation  roads  within  airport 

1 1 .  Service  area  for  rental  cars 

1 2.  Parking  lots  for  rental  cars 

1 3.  Curb  space  for  unloading 

14.  Curb  space  for  loading 

15.  Waiting  areas  at  gates 

16.  Stations  for  intra-airport  transportation  system 

17.  Aircraft  loading  system 

1 8.  Concessions 

1 9.  Rental  car  desk 

20.  Runway  capacity 

21.  Number  of  gates 

22.  Passenger  information 

23.  Cargo  transfer 

24.  Air  traffic  control  system 

25.  Refuse  removal 

26.  Flight  operations  and  crew  facilities 

27.  Aircraft  service  on  the  apron 
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P.  The  problem  can  be  simplified,  for  example,  if  the  function  F  “factors”  into  two  parts;  that 
is,  if  there  are  two  functions  Fa  and  F^  and  if  the  Dj’s  can  be  split  into  two  groups  such  that 

F(Dj,....,D27)»Fa(A)  +  Fb(B) 

where  A  and  B  are  groups  of  Dj’s  such  that  no  D{  is  common  to  both  A  and  B.  A  and  B 
represent  subproblems  that  can  be  solved  separately.  The  general  goal  is  to  break  the  function 
into  as  many  “factors”  as  possible  such  that  there  is  no,  or  very  little,  interaction  between 
factors. 

The  next  step  is  to  determine  the  interaction  between  all  pairs  of  control  variables,  Dj 
and  Dj,  for  example.  Does  the  behavior  of  Dj  with  respect  to  the  performance  function  F 
depend  on  Dj?  Let  R(j  n  be  the  answer  where  R^jjj  may  take  on  one  of  four  values  as  follows: 

0  =  no  obvious  dependency 

1  =  weak  dependency 

2  ~  moderate  dependency 

3  =  strong  dependency. 

Based  on  White’s  subjective  judgment,  values  for  R(jj)  were  generated  and  appear  in  Fig.  17.  It 
is  assumed  that  R^j  jj  =  R(j  j).  Note  that  the  ordering  of  the  items  in  the  matrix  produced  very 
scattered  data.  The  eye  is  not  able  to  identify  any  striking  organizational  structure. 

The  Bond  Energy  Algorithm  was  applied  using  the  original  matrix  as  a  starting  point 
with  the  objective  of  rearranging  the  rows  and  columns  of  the  matrix  to  obtain  a  better  order. 
The  algorithm  tends  to  push  the  larger  numbers  together  into  clumps  and  favors  large  clumps 
over  smaller  ones.  There  is  no  preferential  orientation  of  the  final  clumps;  however,  the 
symmetry  of  the  original  matrix  about  its  diagonal  results  in  a  symmetrical  final  arrangement. 
The  improved  ordering  is  shown  in  Fig.  1 8.  (The  algorithm  applied  to  the  original  matrix 
required  about  2  minutes  of  CDC  1604  computer  time  for  a  number  of  starting  points.) 

After  studying  the  matrix  in  Fig.  18,  it  appeared  that  there  were  eight  clumps  of 
numbers  as  indicated  in  the  figure.  The  clumps  contain  all  of  the  strong  dependencies  (the  3s) 
and  all  but  six  of  the  moderate  dependencies  (the  2s). 

The  interpretation  of  Fig.  18  is  that  clumps  along  the  diagonal  correspond  to  natural 
divisions  of  the  big  problem  into  subtasks.  The  off-diagonai  elements  not  included  in  any 
clump  correspond  to  coordination  links.  Figure  1 9  illustrates  this  interpretation.  The  items  are 


CONTROL  VARIABLES 

I  2  3  4  5  6  7  8  9  10  H  12  13  14  15  16  17  18  19  20  21  22  23  24  2$  26  27 


1  3  3 

2  3  3 

3  0  0 

4  2  3 

3  2  0 

6  0  0 

7  1  I 

8  I  1 

9  0  0 

10  0  0 

'10  0 
12  0  0 

13  2  2 

14  0  0 

15  0  0 

16  0  1 

17  0  0 

18  1  0 

19  0  1 

20  0  0 

21  0  0 

22  3  3 

23  0  0 

24  0  0 

25  0  0 

26  0  0 


2  0  11 
0  0  11 
10  11 
0  0  0  0 
3  0  2  2 
0  3  0  0 
2  0  3  2 

2  0  2  3 
12  2  2 
113  2 
0  0  10 
0  0  0  0 
0  0  I  0 
0  0  10 
0  0  0  0 

3  0  11 
0  10  0 
0  0  0  0 
10  0  0 
0  1  0  0 
3  0  0  0 
10  10 
0  3  0  0 
0  0  0  0 
0  0  0  0 
0  10  0 
0  2  0  0 


0  0  0 
0  0  0 
0  0  0 
0  0  0 

1  I  0 

2  1  0 

2  3  I 
2  2  0 

3  3  0 
3  3  1 
0  1  3 
0  1  3 
0  2  0 
0  2  0 
0  0  0 
0  0  0 
0  0  0 
0  0  0 
0  1  0 
0  0  I 
0  0  0 
0  0  0 
0  0  0 
0  0  0 
0  0  I 
0  0  0 
0  0  0 


0  0  c 
0  0  1 
3  0  1 

1  0  0 
0  0  3 
0  0  0 
1  0  1 
0  0  1 
0  0  0 

2  0  0 
0  0  0 
I  0  0 
0  0  1 
3  0  1 
0  3  1 
1  1  3 
0  I  1 
0  1  1 
I  0  1 
0  0  0 
0  3  3 
0  1  2 
0  0  0 
0  0  0 
0  0  0 
0  0  0 
0  1  0 


0  1  0 
0  0  1 
0  0  1 
1  0  0 
0  0  1 
I  0  0 
0  0  0 
0  0  0 
0  0  0 
0  0  1 
0  0  0 
0  0  I 
0  0  1 
0  0  1 
I  1  0 
1  1  1 
3  0  0 
0  3  0 
0  0  3 
1  0  0 
1  0  0 


0  0  0  0 
0  0  0  0 
0  0  0  0 
0  0  0  0 
0  0  0  0 
0  0  12 
0  0  0  0 
0  0  0  0 
0  0  0  0 
0  0  0  0 
0  10  0 
0  0  0  0 
0  0  0  0 
0  0  0  0 
0  0  0  1 
0  0  0  0 
0  0  0  3 
0  2  0  0 
0  0  0  0 
3  0  10 
10  11 
0  0  0  0 
0  0  0  1 
3  0  2  1 
0  3  0  3 
2  0  3  0 


(  SEE  TABLE  3  EOR  DESCRIPTION  OF  VARIABLES  BY  NUMBER. ) 

FIGURE  17.  Dependency  Matrix  for  Airport  Variables 


CONTROL  VAR.5RUS 

18  25  27  23  6  17  15  21  16  5  8  9  10  7  13  22  1  2  4  3  14  19  12  II  20  24  26 


3  2  0 

2  3r1l 
033! 

0  0  1 

002 

0  0  Ll 

1  0  1 

0  0  1 

1  0  0 

000 

000 

000 

000 

000 

000 

1  0  0 

1  0  0 

000 

000 

000 

000 

000 

000 

0  1  0 

000 

0  0  1 

000 


0  0  1 
0  0  1 
0  0  1 
000 
0  0° 
000 
010 
000 
000 
000 
000 
000 
0  0  1 
000 
000 
000 
000 
000 
0  1  1 
000 
0  1  0 


1  0  t 
000 
1  1  0 
000 
000 
1  1  i 
3  3  1 
3  3  3 
1  3  3 

0  3  3 
0  0  1 
0  0  0 
0  0  0 
0  0  1 
0  0  I 

1  0  d 
0  0  0 
0  0  I 
0  0  0 
0  0  I 
0  0  I 
0  0  I 
0  0  0 
0  0  0 
000 
0  I  0 
0  I  0 


0  0  0  0  0 

0  0  0  0  0 

0  0  0  0  0 

0  0  0  0  0 

0  0  m  1  0 

0  0  0  0  0 

o]o  0  0  0 

3  0  0  0  0 
3  10  0  1 

Tji  i  i  r 

2  3  2  2  2 

1  2  3  3  2 

1  2  3  3  3 

2  2  2  3  3 

6  6  0  LU  1 

10001 
T]  1  0  0  1 
01001 
00000 
1  1  0  0  1 
000  [2]  t 
10010 
00010 
0  0  0  1  1 
00000 
00000 
00000 


0  t  1 
000 
000 
000 
000 
000 
0  1  0 
000 
1  mo 
0  its 
0  0  1 
000 
T]o  0 

1  1  1 
3  2  2 

2  3  3 
2  3  3 
2  3  3 
I  0  2 
0  0  0 
0  0  0 
l  1  0 
1  0  0 
0  0  0 
0  0  0 
0  0  0 
0  0  0 


0  0  0 
0  0  0 
0  0  0 
0  0  0 
0  0  0 
0  I  0 
0  0  0 
0  0  0 
I  0  1 
0  0  I 
I  0  I 
0  0  0 
0  0  0 

1  0  1 

2  1  I  0 

3  0  0 
3  2  0 

3  3  _0_ 
3  [3  13 
0  3  3 
0  1  3 
1  0  I 
0  0  0 
0  0  0 
0  0  0 
0  0  0 
0  0  0 


0  0  0  0  0 
0  10  0  0 
0  0  0  10 
0  0  0  0  0 
0  0  10  1 
0  0  10  0 
0  0  0  0  0 
0  0  [2]  1  I 
0  0  0  0  0 
0  0  0  0  0 
0  0  0  0  0 
0  0  0  0  0 
I  I  0  0  0 
0  10  0  0 
1  0  0  0  0 
0  0  0  0  0 
0  0  0  0  0 
0  0  0  0  0 
0  0  0  0  0 
0  0  0  0  0 
I  0  0  0  0 
I  0] 000 
3  3  0  0  0 
3  3  10  0 
6  1*3  3  I 
0  0  3  3  2 
0  0  12 


FIGURE  18.  Reordered  Dependency  Matrix  for  Maximum  Clumpiness 
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CONCESSIONS - 

REFUSE  REMOVAL - 

AIRCRAFT  SERVICE  ON  APRON - 

CARGO  TRANSFER - 

CARGO  TERMINAL - 

AIRCRAFT  LOADING  SYSTEM - 

WAITING  AREAS  AT  GATES - 

NUMBER  OF  GATES - 

STATIONS  FOR  INTRA-AIRPORT - 

TRANSPORTATION 
INTRA-AIRPORT  TRANSPORTATION 
SYSTEM 

REMOTE  PARKING  LOTS - 

MAIN  ACCESS  ROADS - — 

CIRCULATION  ROADS - 

CLOSE-IN  PARKING  LOTS - 

CURB  SPACE  FOR  UNLOADING - 

PASSENGER  INFORMATION - 

PASSENGER  CHECK-IN - 

BAGGAGE  CHECK-IN - 

BAGGAGE  MOVING  SYSTEM - 

BAGGAGE  CLAIM - 

CURB  SPACE  FOR  LOADING - 

RENTAL  CAR  DESK - 

PARKING  LOTS  FOR  RENTAL  CARS- 
SERVICE  AREA  FOR  RENTAL  CARS  — 

RUNWAY  CAPACITY - 

AIR  TRAFFIC  CONTROL  SYSTEM - 

FLIGHT  OPERATIONS  AND  CREW  — 
FACILITIES 
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FIGURE  19.  Suggested  Subproblems  with  Coordination  Links 


listed  along  the  left  of  Fig.  1 9  in  the  order  found  on  the  “clumped”  matrix  (Fig.  1 8).  As  a  first 
approximation  (shown  below),  the  performance  function  can  be  split  into  eight  factors 
corresponding  to  the  eight  subproblems  shown  in  Fig.  1 9. 

p  85  Fa  (D18>  D25’  D27> 

+  Fb(D27>  d23»  °6’ °17) 

+  Fc(D15,  d21,d16,d5) 

+  Fd  (D5,  Dg,  D9,  D10,  D?) 

+  Fe(D13,’ D22’ Dl>  °2>  D4^ 

+  Ff(D4,  D3,  D14) 

+  Fg(D12-Dll> 

+  Fh  >D20)  D24,  D26) 

Except  for  D27  (aircraft  service  on  apron)  and  D$  (intra-airport  transportation  system),  the 
eight  components  in  the  above  approximation  form  independent  subproblems.  The  six  coordi¬ 
nation  liit'-s  shown  in  Fig.  1 9  could  form  the  basis  of  six  “correction  factors”  which  would 
improve  the  approximation.  The  correction  factors  would  be  of  the  form  Aa  (Dg,  D9),  Ag 
(Dio,  D13>>  Ac(D10’D14)>  Ad(D16»D22^  Ae  (D5»  Di)>  and  Af  (D2q,  f>2l)- 

4.  Ordering  of  Error  Matrices  in  the  Analysis  of  Perception  of  Consonants 

In  this  example,  the  Bond  Energy  Algorithm  is  used  to  reorder  an  error  matrix 
obtained  from  an  experiment  testing  the  perception  of  consonants.  The  matrix  under  consider¬ 
ation  is  a  square  matrix  with  the  29  consonants  lying  on  the  vertical  and  horizontal  axes.  These 
data  were  taken  from  an  article  by  Ahmed  and  Agrawal  (Ref.  29)  in  the  Journal  of  the 
Acoustical  Society.  In  the  experiment  each  consonant  was  enunciated  in  the  intial  position  of 
540  nonsense  syllables.  The  aj  element  of  the  matrix  contained  the  number  of  times 
consonant  0  was  heard  when  consonant  a  was  spoken.  It  is  clear  than  since  the  correct 
consonants  are  heard  most  often,  the  diagonal  elements  of  the  matrix  will  be  largest.  For  this 
example,  the  square  roots  ,4of  the  elements  were  used  rather  than  the  elements  themselves, 
and  all  the  elements  whose  values  are  less  than  two  are  deleted.1 5  The  error  matrix  was  input 
in  a  random  manner  and  the  best  ordering  of  the  reordered  error  mutrix  is  shown  in  Fig.  20. 
The  square  blocks  lying  along  the  principal  diagonal  of  the  matrix  indicate  that  the  consonants 
have  been  clustered  into  7  groups.  These  clusters  represent  those  groups  of  consonants  that 
were  most  often  confused  with  one  another  during  the  experiment.  The  off-diagonal  non-zero 


14.  In  thin  example  th«  weighting  U  again  used  to  ptewve  teak. 

15.  The  unall  dementi  ate  deleted  to  that  the  pattern!  farmed  by  the  target  dementi  may  be  viiually  Identified  mote 
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entries  represent  consonants  in  one  group  being  mistaken  for  consonants  outside  their  group. 
Note,  for  example,  that  one  cluster  contains  the  consonants  d*1,  d,  .d,  and  b  which  occur 
together  because  they  sound  so  much  alike  and  hence  were  often  mistaken  for  one  another 
during  the  experiment. 

This  example  again  illustrates  how  this  method  of  direct  analysis  can  be  a  significant 
aid  in  determining  inherent  group  structure  contained  in  data  matrices. 


Tf  d-dj' tfh  /  i  kh  bh  gh  ph  t*1  .f11  .dh  d*1  d  .d  b  w  r  t  p  g  k  n  i  m  n  I 
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FIGURE  20.  Reordered  ( HINDI )  Consonant  Error  Matrix 


5.  Inter-City  Distances 


In  this  example,  the  Bond  Energy  Algorithm  is  applied  to  a  geographical  problem  to 
determine  if  it  can  satisfactorily  cluster  or  clump  together  neighboring  cities  when  a  large 
number  of  inter-city  distances  are  given  as  input  data.  Since  the  algorithm  clusters  large 
elements  of  a  matrix,  it  was  decided  that  the  square  root16  of  inverse  distance  would  be  used 
for  the  matrix  elements.  In  particular,  the  elements  of  the  matrix  a^^j  are  given  by  the 
expression,  (with  k  =  2) 


-1/2  _  /ioo 

/d" 

i  £  j,  and 

V  d 

a1/2jj=  20 

all  i  =  j 

where  dy  is  the  distance  between  city  i  and  city  j  in  hundreds  of  miles  and  the  resulting  matrix 
elements  are  rounded  to  the  nearest  integer.  This  input  matrix  is  given  in  Fig.  21 ,  where,  for 
visual  clarity,  all  the  elements  with  values  less  than  7  have  been  deleted.  It  should  be 
remembered  from  the  definition  that  the  larger  the  matrix  element,  the  closer  are  the  two 
cities  i  and  j. 


The  matrix  given  in  Fig.  22  is  the  reordered  inverse  distance  matrix  following  operation 
by  the  algorithm  on  the  input  matrix.  Elements  whose  values  are  less  than  7  have  again  been 
deleted. 


A  number  of  clusters  may  easily  be  determined  by  identifying  the  square  blocks  of 
data  that  occur  along  the  main  diagonal  of  Fig.  22.  The  first  two  cities,  Helena,  Montana,  and 
Bismark,  North  Dakota,  are  well  isolated  and  constitute  two  separate  clusters  themselves.  The 
next  two  cities,  Denver,  Colorado,  and  Cheyenne,  Wyoming,  are  quite  close  and  constitute  a 
cluster.  The  next  three  cities,  Des  Moines,  Iowa;  Dubuque,  Iowa;  and  Chicago,  Illinois,  are 
contained  in  the  next  cluster,  and  so  forth.  The  one  anomaly  that  does  exist  is  the  occurrence 
of  the  rectangular  off-diagonal  block  of  7s.  This  indicates  that  although  Chicago,  Detroit,  and 
Ft.  Wayne  are  geographically  near  each  other  and  are  therefore  in  the  same  cluster,  that 
Detroit  and  Ft.  Wayne  are  also  near  some  cities  in  another  cluster,  i.e.,  Cleveland,  Akron, 
Columbus,  and  Cincinnati. 

All  these  clusters  may  be  verified  geographically  by  referring  to  the  map  of  the  United 
States  given  in  Fig.  23.  The  cities  under  consideration  are  denoted  by  darkened  squares  and  the 
clusters  are  shown  by  the  cities  contained  within  each  closed  line. 


16,  llw  square  root  Mfas  used  to  pfeterve  wale  and  keep  the  matrix  elements  teat  than  ot  equal  to  20. 
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The  conclusion  that  can  be  drawn  from  this  example  is  that  the  Bond  Energy 
Algorithm  can  indeed  rearrange  data  geographically  when  it  is  presented  in  another  order 
(alphabetically). 


Preceding  page  blank 


47 


II.  THE  MOMENT  ORDERING  ALGORITHM 


A.  INTRODUCTION 


The  purpose  of  the  Moment  Ordering  Algorithm  is  to  use  the  information  contained  in 
an  array  of  data  to  find  a  one-dimensional  ordering  of  the  row  (and  columns)  of  the  array. 
This  one-dimensional  ordering  will  represent  the  ranking  of  the  rows  (and  columns)  under  the 
relationship  which  the  algorithm  finds  to  be  the  most  important  in  analyzing  the  array.  The 
algorithm  therefore  provides  a  method  of  extracting,  from  the  complex  interrelationships 
which  may  be  expressed  in  the  array,  a  single  important  relationship,  and  of  organizing  the 
rows  and  columns  according  to  this  relationship.  For  example,  one  of  the  problems  discussed 
below  involves  an  array  describing  the  voting  pattern  of  Senators.  The  algorithm  in  this  case 
takes  the  array,  originally  in  the  arbitrary  form  of  an  alphabetical  listing  of  Senators  and  a 
chronological  listing  of  votes,  and  produces  an  ordering  of  the  Senators,  and  of  the  bills  voted 
upon,  based  solely  upon  the  original  array,  which  represents;!  li'  -ral/conservative  ordering.  A 
second  example  involves  an  array  consisting  of  archeological  sites  as  the  columns,  and  of 
pottery  types  rs  the  rows,  with  the  entries  being  the  concentration  of  a  pottery  type  in  a  site. 
The  algorithm  in  this  case  privides  a  reordering  which  puts  the  pottery  types,  and  the  sites,  in  a 
choronological  order,  based  upon  the  fact  that  the  most  important  factor  in  determining  the 
types  of  pottery  found  at  these  sites  was  the  age  of  the  site. 

B.  THE  ALGORITHM 


I .  Motivation 


The  definition  of  the  algorithm  is  based  upon  the  fact  that  if  two  rows  are  similar  to 
each  other,  their  mean  row  moments  should  be  close  to  each  other  in  value.  The  mean  row 


moment  Xj  of  row  i  is  defined  as 
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where  ay  is  the  ij^1  entry  in  the  array.  This  is  merely  another  way  of  stating  that  rows  are 
similar  if  their  large  entries  occur  in  the  same  columns,  or  in  columns  close  to  each  other. 
Similarly,  if  two  columns  are  closely  related,  their  mean  column  moments,  defined  as 


M 

“Hi 

l  =  1 


yi  M 


1  =  I 


should  be  close  to  each  other  in  value. 


Based  upon  these  observations,  then,  it  is  desirable  to  arrange  an  array  so  that  its  rows 
are  in  order  of  the  values  of  their  row  moments,  while  at  the  same  time  its  columns  are  in 
order  of  the  values  of  the  column  moments.  This  state  will  correspond  to  a  one-dimensional 
ranking  of  both  the  rows  and  the  columns  according  to  the  same  underlying  variable.  The 
algorithm  provides  a  method  of  finding  such  states,  and  hence  of  ordering  arrays  of  data. 

2.  Definition 

The  algorithm,  beginning  with  an  arbitrary  arrangement  of  an  array,  proceeds  in  the 
following  way  to  find  a  state  with  the  property  described  above,  of  having  both  the  rows  and 
the  columns  of  the  array  arranged  in  order  of  their  moments: 

1 .  The  row  moments  are  calculated  for  the  original  arrangement  of  the  array,  and 
the  rows  are  reordered  to  put  them  in  order  of  their  moments. 

2.  The  column  moments  are  calculated,  and  the  columns  reordered  according  to 
their  moments. 

3.  Because  the  reordering  of  the  columns  changes  the  values  of  the  row  moments, 
the  rows  will  no  longer  necessarily  be  in  order  of  their  row  moments.  The  row 
moments  are  therefore  recalculated  for  the  new  arrangement  of  the  columns, 
and  the  rows  reordered  according  to  these  new  moments. 

4.  The  procedure  is  continued,  alternately  reordering  the  rows  and  columns,  until 
a  state  is  found  in  which  both  are  simultaneously  in  order  of  their  moments. 
This  state,  then,  is  the  desired  ordering  of  the  rows  and  columns,  and  is  a 
solution  of  the  algorithm. 


SO 


The  algorithm  is  therefore  entirely  an  iterative  procedure.  The  progress  of  the  algo¬ 
rithm  toward  convergence,  however,  is  marked  by  an  increasing  concentration  of  the  larger 
elements  on  or  near  the  main  diagonal.1 

The  progress  of  the  algorithm  is  illustrated,  for  a  4x4  array,  in  Fig.  24.  The  initial  state 
of  the  array  is  a;  the  values  of  the  row  moments  for  that  array  arrangement  are  also  shown. 
The  algorithm  then  proceeds  through  states  b,  c,  and  d,  by  reordering  the  rows  and  columns 
alternately.  When  state  e  is  reached,  it  is  found  that  the  rows  are  already  in  the  proper  order 
and  do  not  need  to  be  reordered.  This  marks  that  state  as  a  solution. 

The  concentration  of  the  larger  elements  on  or  near  the  main  diagonal  in  the  solution  is 
pointed  out  in  Fig.  24  by  circling,  in  the  initial  and  final  states,  the  four  largest  elements.  They 
arc  scattered  in  the  initial  state  but  in  the  solution  three  are  on  the  main  diagonal  and  one  is 
just  off  it. 

The  following  subsections  present  further  details  concerning  the  use  of  the  algorithm. 
Section  C  presents  several  specific  problems  which  have  been  investigated  by  use  of  the 
algorithm,  and  illustrates  the  utility  of  the  orderings  produced  by  the  algorithm. 

3.  Stable  States  and  Multiple  Solutions 

The  algorithm  as  defined  above  takes  an  arbitrary  initial  ordering  of  an  array  and  finds 
a  stable  ^orderin'1.  It  has  been  found,  however,  that  if  different  initial  orderings  of  the  same 
array  are  used,  different  solutions  may  be  found.  For  example.  Fig.  25  shows  two  different 
solutions  whicn  can  be  found  for  a  simple  3x3  array.2 

When  the  algorithm  is  run  many  times  on  larger  arrays,  using  different  starting 
orderings  each  time,  it  has  been  found  that  those  solutions  which  occur  most  frequently 
always  are  amongst  the  most  diagonal  arrangements  of  the  array.3  Conversely,  any  solutions 
which  are  very  nondiagonal  occur  only  rarely. 

This  observation  has  been  used  as  the  basis  of  a  technique  for  obtaining  a  final  ordering 
of  the  rows  and  columns  which  best  utilizes  the  additional  information  found  in  the  multiple 
solutions. 

Appendix  H  discusses  a  measure  of  effectiveness  which  has  been  defined  to  measure  this  progress  toward 
dtagoitatity.  However,  because  unlike  the  Bond  Energy  Algorithm  this  algorithm  was  not  developed  tomaumize  this  quantity, 
the  measure  of  effectiveness  defined  has  been  found  to  be  of  only  marginal  use. 

2.  Appendix  I  describes  an  investigation  which  was  made,  for  a  3x1  array,  of  the  properties  which  lead  io  the  existence 
of  these  multiple  solutions. 

3.  As  measured  by  the  correlation  coefficient  measure  of  effectiveness  defined  in  Appendix  H. 
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FIGURE  24.  Operation  of  the  Algorithm  on  a  Small  Array 
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FIGURE  25.  Illustration  ot  Multiple  Solutions 


The  algorithm  is  run  a  “large”  (25  or  50  has  been  found  satisfactory)  number  of  times, 
each  time  starting  from  a  different  random  ordering  of  the  rows  or  columns,  and  the  order  of 
the  rows  and  columns  found  each  time  is  saved.4  The  average  of  the  position  taken  by  each 
row  (and  column)  in  the  solutions  is  found.  (Solutions  found  more  than  once  are  entered  once 
for  each  time  found  in  obtaining  the  average.)  The  rows’  and  columns’  final  order  is  then 
simply  the  order  of  their  average  positions.  Most  often,  this  order  will  be  the  same  as  the  order 
in  the  most  common  solution;  it  is  always  very  close  to  that  order. 

Despite  the  additional  complication  introduced,  this  technique  is  considered  preferable 
to  merely  taking  the  most  common  solution,  because  in  the  event  that  several  solutions  are 
common,  this  technique  best  takes  into  account  the  alternative  orderings  each  solution 
represents  in  arriving  at  a  consensus  final  ordering. 

4.  Additional  Details 

The  previous  sections  have  discussed  all  of  the  features  of  the  algorithm  which  are 
important  in  pract'ce.  There  are,  however,  two  points  of  theoretical  interest  which  must  be 
mentioned  at  this  poin'  Both  concern  situations  which  can  arise  in  the  process  of  iteration 
carried  out  by  the  algorithm.  Both  occur  so  rarely,  however,  that  in  practice  they  can  usually 
be  ignored. 

a.  Ties.  In  carrying  out  the  algorithm,  two  or  more  rows  or  columns  may  have 
identical  moments.  In  this  case  it  is  necessary  to  resolve  the  tie  to  obtain  an  ordering  so  that 
the  algorithm  can  proceed.  This  is  done  by  trying  all  permutations  of  non-identical  rows  (or 
columns)  and  selecting  that  particular  row  (or  column)  order  which  yields  the  highest  value  of 


4.  Note  that  a  particular  order  and  Us  reverse  are  considered  identical  and  saved  as  the  same  order. 
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the  correlation  coefficient  R.5  If  several  permutations  of  the  rows  (or  columns)  have  the  same 
value  of  R,  the  algorithm  simply  accepts  the  last  order  investigated.  It  should  be  noted  that 
ties,  while  prominent  for  small,  binary  (0-1 ),  arrays,  very  rarely  occur  when  dealing  with  large 
arrays  containing  non-binary  data. 

b.  Cycling.  According  to  the  definition  of  the  algorithm,  the  iterative  procedure  is 
continued  until  a  stable  state  unchanged  by  either  row  or  column  operations  (Fig.  24,  for 
example)  is  found.  In  fact,  however,  it  is  theoretically  possible  that,  instead  of  arriving  at  such 
a  stable  state,  the  algorithm  may  cycle  between  a  small  set  of  states.  Fig.  26  illustrates  the 
phenomenon  for  a  specially  designed  small  array  (in  actual  fact  such  cycling  has  only  been 
found  in  very  much  larger  arrays).6  Once  the  algorithm  arrives  at  the  state  shown  in  Fig.  26, 
which  it  can  reach  from  many  other  states,  it  will  cycle  forever  between  a,  b,  c,  and  d,  in  that 
order.  Such  an  “infinite  loop"  itself  represents  a  final  state  of  the  array.  The  procedure  used 
when  such  cycling  is  detected  therefore  is  to  terminate  the  iterations  and  take  one  of  the 
states  involved  in  the  loop  as  the  solution. 
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FIGURE  26.  Illustration  of  Cycling  Phenomenon 

In  practice,  this  phenomenon  has  been  observed  only  very  rarely,  and  only  in  very  large 
arrays.  Furthermore,  even  when  it  does  occur,  it  has  been  found  that  most  often  the  algorithm 
will  find  normal  stable  solutions  when  operating  upon  the  same  array  from  other  starting 
points.  For  this  reason,  this  cycling,  while  theoretically  quite  objectionable,  has  been  found  to 
be  of  little  operational  difficulty. 

5.  Sec  Appendix  H. 

6.  The  symmetry  and  normalization  inherent  in  the  array  of  Fig.  26  ore  not  necessary  for  the  cycling  to  occur,  but 
wero  built  in  to  simplify  the  array. 
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5.  Computation  Time 


For  an  MxN  matrix,  M  operations  are  required  to  compute  each  column  moment  and 
N  operations  to  compute  each  row  moment.  Therefore,  for  each  iteration,  the  total  number  of 
operations  necessary  to  reorder  the  rows  and  columns  is  2MN.  Finally,  if  it  requires  1  iterations 
for  the  algorithm  to  converge,  the  total  number  of  operations  to  reach  a  solution  for  each 
random  starting  point  is  2IMN.  The  computer  time  required  on  the  CDC  1604  to  solve  a 
particular  29x29  matrix  was  about  24  seconds  for  one  starting  point;  an  80x80  matrix  took  4 
minutes.  Note  that  these  times  are  influenced  by  the  number  of  iterations  required  for 
convergence  as  well  as  by  the  matrix  sizes. 

C.  RESULTS 

This  section  describes  two  problems  investigated  with  the  Moment  Ordering  Algorithm. 
It  demonstrates  that  the  algorithm  can  in  fact  uncover  a  dominant  relationship  from  the  vast 
amount  of  information  in  a  matrix,  and  can  produce  orderings  of  the  rows  and  columns  which 
reflect  this  relationship. 

1 .  U.S.  Senate  Voting  Patterns 

The  algorithm  was  used  to  study  the  relationships  between  the  voting  patterns  of  a 
group  of  U.S.  Senators.  The  hope  was  that,  given  only  the  recorded  positions  of  Senators  on  a 
random  group  of  issues,  the  algorithm  could  generate  a  meaningful  ordering.  The  First  20 
Senators  (alphabetically)  in  1968  were  chosen,  and  their  recorded  positions7  on  12  issues  were 
tabulated  (see  Tables  4  and  5).  The  recorded  position  of  the  President  on  each  issue  was  added 
to  the  table,  and  the  algorithm  was  applied  to  the  resulting  21x12  array.  The  results,  as  shown 
in  Table  6,  showed  an  ordering  from  conservative  Republican  and  Southern  Democrat  at  one 
end  to  liberal  Democrat  on  the  other.  To  be  sure  that  the  strong  ordering  was  not  an  accident, 
the  same  type  of  array  was  constructed  tor  12  different  roli  calls  (but  the  same  Senators),  and 
the  algorithm  was  rerun.  The  correlation  between  the  two  sets  of  results  (sec  Table  6  again) 
indicates  that  the  ordering  found  was  significant.  The  difference  between  the  two  rankings 
does  not,  it  is  emphasized,  reflect  any  inherent  limitation  upon  the  accuracy  of  the  algorithm, 
but  rather  is  a  result  of  the  limited  sizes  of  the  samples  of  votes  used  in  the  analyses.  If  more 
roll  calls  were  added  to  the  arrays,  the  results  would  approach  each  other  more  and  more, 
reflecting  the  enlarged  and  therefore  improved  sampling.  The  algorithm’s  solution  indicates 
that,  as  might  be  expected,  although  a  Senator’s  position  on  any  given  issue  may  not  always  be 


7.  As  taken  from  tables  in  Concessional  Quarterly  Almanac,  Voi.  24, 1968,  pp.  IS-58S. 


55 


SENA  TE  ROLL  CALL  VO TES INCL  UDED  IN  ANAL  YSIS 


Subjects  of  roll  calk  identified  in  fable  5. 
Republicans  in  capital  letters. 


Table  5.  ARRA  YS  USED  IN  SENA  TE  VOTE  PA  TTERN  ANAL  YS/S * 


Vote 

No. 

Roll*1 

Call 

Vote 

Subject  Matter 

Sponsor 

10 

jj 

58 

Amendment  to  open  housing  bill  to  bar  federal  courts  from  impairing  title  to  real 
property  as  recorded  und*r  state  recording  statutes. 

Fnun 

■s 

20 

6’ 

21 

Amendment  to  open  housing  bill  to  punish  anyone  instructing  m  the  use  of  fire 
arms  for  riots,  or  interfering  with  police  during  a  riot. 

Long 

3 

30 

61 

19 

Amendment  to  open  housing  bill  to  provide  a  compromise  bill. 

Dirksen 

4 

40 

19 

5K 

Amendment  to  gold  cover  removal  hill  to  limit  cspansion  of  Federal  Reserve 
notes  in  circulation  ro  4'-'  per  year. 

Allot! 

5 

50 

43 

28 

Amendment  to  Standards  of  Conduct  Resolution  r<»  allow  use  of  political 
contributions  for  certain  office  expenses. 

Yarborough- 

Jasits 

6 

60 

3k 

44 

Amendment  to  excise  tax  extension  bill  to  provide  2(71  surtax  on  people 
trading  with  Communist  nations  which  supply  North  Vietnam 

Mundt- 
ByrdiVa  i 

7 

70 

53 

35 

Amendment  to  excise  tax  extension  bill  to  impose  10*'  income  tax  surcharge 
and  limit  expenditures  to  $1 80  billion. 

Wdlunis- 

Smathvs 

K 

»o 

:* 

30 

Amendment  to  Military  Procurement  Authorization  to  cut  RAI>  funds  trout 

57  <>  to  57  4  billion. 

Hart 

*) 

90 

39 

29 

Amendment  to  Conservation  Fund  bill  to  remove  outer  continental  shelt 
revenues  Irom  turn!  for  1972  and  1973 

Williams 

10 

100 

:o 

S3 

Amendment  »o  Omnibus  Crime  Bill  to  prohibit  interstate  mailorder  vdes 
ot  nflc*s  and  shotguns 

kerineds 

It 

no 

51 

30 

Amendment  to  Omnibus  Cnme  Bill  to  delete  language  dcnvng  Supr»  <  ourt 
jurisdiction  to  review  state  court  judges’  decisions  to  admit  s  >evmness  testimony 
in  evidence 

1  s  dings 

1? 

i:o 

5.1 

44 

Amendment  t«  -  -thnmhux Crime  Bill  to  allocate  2 '3  instead  of  85  •'  ot  funds  m 
block  grants  to  states. 

Brooke 

13 

i:?1- 

25 

35 

Amendment  to  require  cities  as  well  as  states  to  reimburse  NIIK  lor  vosts  ot 
h*»l  losses  insured  bv  MIX’. 

Russell 

14 

1 57s 

42 

:? 

Amendment  to  delete  section  on  retirement  benefits  from  bill  to  extend  term 
ot  office  nl  bankruptcy  relerees 

Carlson 

150 

44 

32 

Motion  to  table  amendment  winch  would  have  provided  552  million  supplc- 
mental  anpropri.it ion  to  Labor  Depuftmcnl  for  summer  jobs. 

Holland 

lh 

i  60 

16 

61 

Amendments  to  Military  Construction  Authorization  to  ml  Navy  ami  An 

Force  funds  by  10  • 

Clark 

17 

i7o 

34 

tx 

Aiuemftnenl  io  luvenue  delinquency  bill  to  allocate  all  funds  j\  block  grants  to 

stales. 

Murphy 

1* 

U40 

34 

52 

Amendment  l«i  f  vderal  Agency  Authorization  to  cut  NASA  HAD  lunds  an 
additional  5  t(M)  nuliion 

Williams 

*•» 

1% 

30 

40 

Amendment  to  Agricultural  Act  to  limit  to  575.000  payments  in  one  prodwr 
for  participation  in  certain  agricultural  programs. 

Williams 

20 

200 

46 

45 

Amendment  to  Interest  Rales  Hill  to  strike  out  language  authorizing  Federal 
Reserve  hanks  to  purchase  obligations  directly  from  federal  agent  lev 

Hennrtt 

:i 

:nl 

51 

S  1 

Amendment  to  strike  mil  language  added  by  House  which  limited  expenditures 
of  State.  Justice  and  Commerce  to  SI  98  billion 

Committee 

:: 

221" 

46 

28 

Foreign  Aid  Authorization  Hill 

— 

M 

230 

23 

35 

Amendment  to  Renegotiation  Act  lu  exempt  Renegotiation  Hoard  from 
employee  limitations- 

Proxmire 

24 

M0 

31 

53 

Amendment  to(iun  Control  Act  to  add  a  registration  provision. 

Brooke 

d  Information  is  taken  from  the  (imgreiuonat  Quarter!'  AlmantH'.  iVol.  '4.  196X1.  pp  I.V5KS. 

**  fhe  first  1 2  roll  calls  above  are  included  in  the  first  array,  the  second  1 2  in  the  second  array 

1  When  a  roil  call  selected  was  too  iine-sided  to  convey  ogmfkant  information,  a  roll  call  close  in  time  to  H  was  substituted 
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CON 


Table  6.  RESUL TS  OF  VOTE  PA  TTERN  ANAL  YSIS • 


Order 

Array  lb 

Array  2b 

1 

Burdick 

Burdick 

i  2 

Clark 

Bayh 

3 

Church 

Clark 

4 

Brewster 

Brewster 

5 

Bible 

Church 

6 

CASE 

Bartlett 

j  7 

Bayh 

Byrd  (W.Va.) 

j  8 

Anderson 

Anderson 

3  9 

Bartlett 

CASE 

1  10 

President 

Cannon 

J 

:  n 

BROOKE 

COOPER 

l  12 

Cannon 

BROOKE 

a  13 

o 

AIKEN 

Bible 

5  >4 

Byrd  (W.Va.) 

President 

15 

ALLOTT 

BOGGS 

16 

BAKER 

AIKEN 

17 

BOGGS 

ALLOTT 

18 

COOPER 

Byrd  (Va.) 

19 

BENNETT 

BAKER 

,  20 

CARLSON 

CARLSON 

BENNETT 

21 

Byrd  (Va.) 

predictable,  overall  voting  patterns  based  upon  ideology  are  strongly  evident,  and  Senators  can 
be  placed  reasonably  well  on  a  liberal-conservative  spectrum.  More  important,  for  our 
purposes,  it  indicates  that  when  a  meaningful  ordering  is  inherent  in  a  set  of  data,  the 
algorithm  will  find  that  ordering. 

2.  Chronological  Ordering  in  Archaeology 

The  algorithm  was  used  to  attempt  to  order  a  series  of  archaeological  deposits.  The 
basic  data  available  is  the  distribution  of  various  types  of  pottery  (eight,  in  this  case)  among 
various  deposits  of  archaeological  interest  (also  eight,  in  this  case).  Robinson  (Ref.  30),  upon 
whose  work  this  example  is  based,  hypothesized  that  it  should  be  possible  to  arrange  these 
sites  into  a  proper  chronological  order  by  assuming  that  pottery  types  come  into  and  go  out  of 
general  use  in  a  regular  manner  over  time,  and  that,  therefore,  deposits  similar  to  each  other  in 
the  amounts  of  various  types  of  pottery  will  be  close  to  each  other  in  time  as  well.  Thus,  if  a 
satisfactory  one-dimensional  arrangement  of  the  pottery  deposits  can  be  found,  on  the  basis  of 
a  pottery-type  percentage  array,  the  sites  should  be  chronologically  ordered.  This  was 
therefore  used  as  a  test  of  the  Moment  Ordering  Algorithm. 

The  raw  data  matrix  presented  by  Robinson  in  Ref.  30  is  shown  in  Table  7.  If  the 
algorithm  is  performed  on  this  array,  the  solution  found  is  3A,  2A.  3B.  I  A,  3C,  2B,  IB.  20. 
which  is  very  close  to  that  presented  by  Robinson,  and  which  satisfies  the  tests  he  carries  out 
on  his  candidate  solution. 


Table  7.  RA  W  POTTER  Y  PLRCENTA  GES 


Pottery 

Deposit 

Type 

2A 

28 

2C 

1A 

IB 

3A 

3B 

3C 

l 

24.0 

1.4 

0.2 

11.3 

0.3 

29.6 

54.3 

0 

2 

66.8 

0.9 

0 

0 

0 

0 

3.5 

0 

3 

1.3 

0 

0.2 

3.8 

0.2 

14.1 

14.0 

6.6 

4 

0 

0 

0 

1.3 

0.2 

0 

1.8 

3.3 

5 

0 

0 

0 

3.3 

0.5 

0 

5.3 

5.5 

6 

4.0 

0 

0 

24.9 

1.4 

7.0 

7.0 

27.5 

7 

0 

97.7 

99.3 

52.6 

97.4 

0 

12.3 

57.1 

8 

3.9 

0 

0.3 

2.8 

0 

49.3 

1.8 

0 
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Robinson,  however,  introduces  an  “agreement  coefficient”  between  two  pottery  types, 
defined  arbitrarily  as: 

N 

ay  =  200  -  £  I  Pjk  -  Pjjc  l  .  where  there  are  N  sites, 
k  =  1 

and  and  Py,  are  the  percentages  of  types  i  and  j  in  site  k.  Therefore,  «  200  constitutes 
total  agreement  between  the  composition  of  two  sites,  Ujj  =  0,  total  disagreement.  Robinson’s 
resulting  array  is  presented  in  Table  8.  Robinson  then  attempts  to  carry  out  a  “rearrangement” 
of  this  array  to  drive  large  numbers  toward  the  diagonal;  he  describes  a  semi-systematic  manual 
method  of  doing  so  and  presents  the  resulting  order  as  his  solution.  The  Moment  Ordering 
Algorithm  was  run  on  Table  8  and  found  exactly  the  same  order  as  Robinson’s  method  -2A. 
3 A,  3B,  1  A,  3C,  IB,  2B,  2C.  The  reordered  matrix  of  agreement  coefficients  is  shown  in  Table 
9,  where  it  is  apparent  that  the  larger  matrix  elements  have  accumulated  around  the  main 
diagonal  of  the  array.  The  advantage  obtained  in  using  the  algorithm,  of  course,  lies  in  the  fact 
that  it  is  an  automatic,  systematic  approach  and  does  not  require  personal  judgments  to  be 
made,  as  Robinson’s  method  did.  The  fact  that  it  reproduces  Robinson’s  chronological 
ordering  reinforces  the  belief  that  the  algorithm  is  suitable  for  just  this  problem -ordering 
entities  in  one  dimension  based  on  their  interrelationships. 


Table  8.  AGREEMENT  COEFFICIENTS 


Pottery 

Pottery 

Deposit 

Deposit 

2A 

2B 

2C 

1A 

IB 

2A 

3B 

3C 

2A 

200 

5 

1 

39 

4 

66 

69 

11 

2B 

5 

200 

196 

108 

195 

3 

29 

114 

2C 

j 

196 

200 

107 

196 

1 

26 

115 

1A 

39 

108 

107 

200 

110 

50 

82 

172 

IB 

4 

195 

1% 

110 

200 

4 

30 

119 

3A 

66 

3 

1 

50 

4 

200 

101 

27 

3B 

69 

29 

26 

82 

30 

101 

200 

66 

3C 

11 

114 

115 

172 

119 

27 

66 

200 
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III.  THE  MOMENT  COMPRESSION  ALGORITHM 


A.  INTRODUCTION 

The  Moment  Compression  Algorithm  discussed  in  this  chapter  is  based  upon  the  key 
observation  that  the  distinguishing  feature  of  a  matrix  in  perfect  block  form,  (see  sketch)  when 
contrasted  with  the  same  matrix  after  row  or  column  permutations,  is  that  the  moment  of 
inertia  of  each  row  and  column  about  its  mean  is  minimized:  any  row  or  column  permutation 
of  a  matrix  in  perfect  diagonal  block  form  will  “expand”  a  block  and  make  it  less  dense, 
thereby  increasing  the  matrix’s  summed  moments  of  inertia. 


Consequently  a  procedure  which  minimizes,  by  row  and  column  permutations,  the 
sums  of  the  row  and  column  mean  square  moments  about  their  means  will  drive  the  matrix 
into  perfect  block  form  if  this  is  possible.1  ,J  if  this  is  not  possible,  the  procedure  will  still  tend 
to  produce  a  pleasing  pattern  because  it  tries  to  create  dense  blocks.  This  reasoning  led  to  the 
development  of  the  Moment  Compression  Algorithm. 

Although  Moment  Compression  has  b«*en  superseded  by  Bond  Encife/  both  as  a 
theoretical  ME  and  as  a  computational  procedure,  this  material  is  being  presented  both  to 
indicate  an  approach  which  was  explored  and  found  impractical,  and  to  show  a  logical 
stepping-stone  in  the  development  of  the  Bond  Energy  Algorithm.  Moment  Compression  was 
historically  important  for  four  reasons: 


will  b«  considered  equally  good.  Bui  one  would  be  indifferent  to  such  ambiguity  as  long  as  Ibc  variables  have  been  factored 
correctly. 

2.  This  assertion  is  proved  in  Appendix  C. 

Preceditig  page  blank 
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( 1 )  It  was  our  first  attempt1 * 3  to  describe  the  appeal  of  a  pattern  in  terms  of 
a  quantitative  ME,  the  sums  of  the  moments  of  inertia.  This  was 
motivated  by  a  desire  to  produce  dense  blocks  of  numbers. 

(2)  It  was  our  first  attempt  to  devise  an  algorithm  based  on  ME- 
optimization.  This  was  in  contrast  to  heuristic  algorithms,  such  as 
moment  ordering  and  some  similarity  matrix  approaches,  where  it  was 
not  clear  what  each  step  in  the  algorithm  was  trying  to  accomplish.  In 
particular,  rigorous  optimization  of  the  ME  would  avoid  the  problems 
of  cycling  and  non-uniqueness4  experienced  in  the  Moment  Ordering 
Algorithm. 

(3)  It  was  our  first  attempt  to  devise  algorithms  which  find  near-optimal, 
rather  than  optimal,  solutions  for  the  ME.  The  major  pitfall  encoun¬ 
tered  in  the  Moment  Compression  case,  but  not  in  the  Bond  Energy 
case,  was  that  the  approximate  algorithm  was  slow5  and  poor.6 

(4)  It  used  an  ME  which  decomposed  into  two  parts,  one  (sum  of  the  row 
moments)  dependent  only  on  column  permutations  and  the  other  (sum 
of  the  column  moments)  dependent  only  on  row  permutations.  Conse¬ 
quently  optimization  of  the  ME  could  be  achieved  in  exactly  two 
passes,  one  finding  the  optimal  column  permutation,  the  other  finding 
the  optimal  row  permutation.  These  two  passes  are  carried  out  com¬ 
pletely  independently  of  each  other,  and  in  particular,  it  is  not 
necessary  to  alternate  between  row  and  column  permutations,  as  in  the 
Moment  Ordering  Algorithm.  This  decomposition  of  the  ME  into  two 
parts  was  an  attractive  feature  later  used  in  the  Bond  Energy  ME 
(row^bonds  and  column-bonds  being  optimized  separately). 


1.  D».  Gould  had  earlier  suggested  use  of  the  matrix  correlation  coefficient  as  a  guide  to  the  performance  of  the 

Moment  Ordering  Algorithm,  but  there  was  no  particular  pattern  that  one  hoped  to  drive  the  matrix  into. 

4.  Cycling  can  never  occur  in  an  algorithm  which  iteratively  optimizes  an  ME,  for  the  ME  is  monotone  from  one 
iteration  to  the  next.  There  would  still  be  non-uniqueness  if  a  few  permutations  achieved  the  global  optimum;  this  could  be 
expected  only  in  degenerate  cases,  and  normally  would  not  occur.  Permutations  achieving  local  (rather  than  global)  optima  of 

the  MF.  could  be  discarded  on  the  basis  of  their  inferior  MEs,  i n  that  many  fewer  “stable"  solutions  could  be  expected  than  in 
the  Moment  Ordering  Algorithm. 

5.  At  least  a  factor  of  three  slower  than  the  Bond  Energy  Algorithm,  and  therefore  impractical  for  problems  larger 
than  about  25x2$. 

6.  While  the  algorithm  is  always  successful  at  putting  a  matrix  into  near  block  diagonal  form,  if  t)  Is  it  possible,  it  had 
two  major  weaknesses  of  (1)  sensitivity  of  the  result  to  the  starting  point,  and  (2)  an  inability  to  handle  the  checkerboard  case, 
shown  in  Fig.  6. 
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B.  MEASURE  OF  EFFECTIVENESS  FOR  MOMENT  COMPRESSION 


As  stated  above,  the  measure  of  effectiveness  for  Moment  Compression  is  the  sum  of 
the  mean-square7  column  moments  and  mean-square  row  moments.  For  any  NxM  non¬ 
negative  matrix  (by),  the  ME  is 

N  M 

ME(b)  =  £  ri  +  £  cj 
i=l  j=l 

where  rj  is  the  row  moment  for  the  i**1  row: 


r  ^  1 

M  I  k?, 

ri,b)  -  n  *  E  »u  !•— - 

j  -  1  M 

L  n  =  1 

M  ” 

£>.  bini 
m  =  1 


and  C:  is  the  column  moment  for  the  jth  column: 


Let  A  =  j^ayj  ie  the  original  NxM  non-negative  matrix  and  let  ^byj  =  [ai  »r(j)J 
denote  the  matrix  whose  j**1  column  is  the  column  of  A,  where  tr  =  J  t ( 1 ),  w(2), . , . ,  ir(V 

denotes  a  permutation  of  1 1,  2 . M  }  .  The  problem  of  finding  the  best  column  permutation 

of  A  is  given  by 

?.  White  any  even  moment  can  be  uted,  tbc  tecond  moment  it  the  umpteu. 


min 


t 

min 

t 


N 

E  1(b) 

i=  1 

M 

jk?i  9flw(i)»00 


where 


N 


Qjkrs  “  ,Si 


1  =  i 


airJ2  Sjk 


W; 


airais*k 


w. 

1 


1<  jkrs  <  M 


M 


and  Wj  -  £  3|j  denotes  the  row  sum  for  the  row. 


j=  1 


Finding  the  best  row  permutation  leads  to  a  problem  completely  analogous  to  that  of  finding 
the  column  permutation. 


The  above  problem  involves  a  minimization  over  all  M!  possible  permutations.  It  is 
called  a  quadratic  assignment  problem  because  of  the  double  appearance  of  it  in  the  minimand. 
The  problem  of  ME  optimization  is  consequently  equivalent  to  solving  two  quadratic  assign* 
mcnt  problems.  (Appendix  A  demonstrates  that  the  same  holds  true  for  the  Bond  Energy  ME.) 

As  discussed  in  Appendix  A,  exact  algorithms  for  solving  quadratic  assignment  prob¬ 
lems  are  too  time  consuming  to  be  practical  for  M  larger  than  1 5  or  20.  Consequently,  an 
approximate  algorithm  was  employed  to  find  a  near-optimal  solution.  The  approximate 
algorithm  is  a  gradient  search  in  M “-dimensional  space,  and  is  described  in  the  next  section. 


C.  GRADIENT  ALGORITHM  FOR  APPROXIMATE  ME  OPTIMIZATION 


The  minimization  problem  posed  above  can  be  rewritten  as 


min  Z(X) 
XePM 


where  PM  denotes  the  set  of  alt  M!  possible  MxM  permutation  matrices  ( i.e.,  all  matrices  of  the 
torm  Xjj  =  and  where 


6* 


M 

Z<X>  =  E  Qjkrsxjrxks  =  <b>x>  +  (X.C.X) 

jkrs  =1 


M  M 

(B’X)=.,E.  BjkXjk  =  E 

jk  =  1 


jk  -  1 
M 

■I 

jkrs=  1 


v  !il£ 

i  =  1  wi 


M  M 

(X.C.X)  =  £  Xjkxrsc,krs  =  £ 

'  jkrs  =  1 


Xjk 


N  aikaisir 
“i=l 


XjkXrs 


Note  that  C  is  negative  semi-definite: 

(y  C.y)  ^  o  for  any  MxM  matrix  y. 

The  gradient  search  was  motivated  by  a  paper  by  Dem’ianov  (Ref.  23)  Exploiting  the 
quadratic  dependence  of  Z  and  the  negative-definiteness  of  C,  one  writes 

Z(X)  =  Z(X°)  +  (X-X°.  grad  Z(X0))  +  (X-X°,C,X-X°) 

where  the  last  term  is  non-positive  and  where 

M 

grad  Z(X°)jk  =  Bjk  +  J  £  Cjkls  X°rs 

rs  =  l 

Consequently  if  X  is  chosen  to  minimize  (X,  graa  Z(X°)).  one  finds  Z  (X)  <  Z  (X°).  with 
equality  usually  implying  that  X°  is  local  minimum  of  Z.8  The  following  gradient  algorithm  is 
the  result: 

Step  1.  Select  an  initial  permutation  matrix  X°'d. 

Step  2.  Compute  grad  Z(X°*d). 

Step  3.  Solve  min  (X,  grad  Z(X°*d))  for  the  minimizing  permutation  matrix 

XePM 

X»ew 

Step  4.  If  Xncw  =£  X°*d,  set  Xob*  B  Xncw  and  return  to  Step  2;  if  Xnew  » 

Xold,  stop. 

The  algorithm  converges  to  a  permutation  matrix  which  generally  is  a  local  minimum  of  Z(X).8 
The  time  consuming  portion  of  the  algorithm  is  Step  3.  The  minimization  in  Step  3  is 


8.  The  basic  property  in  that  if  (X-X0,  yad  2(X°))>Q  for  at]  permutation  matrices  X,  and  if  (X-X°,  C,X-X°)  =  0  for 
aft  X  for  which  there  is  strict  equality,  then  X°  is  a  local  minimum  for  Z(-),  where  the  domain  of  2  is  now  extended  to  the  set 
ot  all  doubly  stochastic  matrices. 
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[M 

T  grad  Z(X°) . 
j=l  > 

where  n  ranges  over  all  permutations  of  {l  ,2, . . . ,  M }  . 

This  class  of  problems  is  known  as  linear  assignment  problems  and  is  most  readily  solved  by 
the  so-called  Hungarian  method  (Ref.  32).  Unfortunately,  the  labor  for  the  Hungarian  method 
is  proportional  to  or  and  since  several  linear  assignment  problems  must  be  solved,  the 
computation  time  for  this  gradient  algorithm  turns  out  to  be  excessive  for  large  M. 

D.  COMPUTATIONAL  RESULTS 

The  gradient  algorithm  described  above  was  coded  in  order  to  provide  near-optimal 
solutions  to  the  Moment  Compression  problem.  The  gradient  algorithm  is  used  twice;  once  to 
minimize  the  sum  of  the  row  moments  and  again  for  the  column  moments.  The  major 
computational  effort  goes  into  solutions  of  successive  linear  assignment  problems. 

The  primary  advantages  of  the  gradient  algorithm  are  its  simplicity  and  (as  the 
following  two  examples  illustrate)  its  excellent  capability  for  putting  a  matrix  into  near 
block-diagonal  form  when  this  is  possible.  The  primary  disadvantage  is  the  large  computer  time 
(a  factor  of  three  greater  than  for  the  Bond  Energy  Algorithm),  rendering  the  method 
impractical  for  matrices  larger  than  ..'out  25x25. 

The  excessive  computational  e  fort  arises  from  two  sources.  One  is  the  need  to  solve 
successive  linear  assignment  problems,  each  of  which  is  time  consuming.  The  second  is  the 
existence  of  several  local  minima  for  Z(X),  with  the  consequence  that  the  final  data  ordering  is 
somewhat  sensitive  to  the  initial  data  ordering.  (The  Moment  Ordering  Algorithm  has  similar 
properties.)  It  therefore  is  necessary  to  start  the  algorithm  at  several  randomly-selected  initial 
permutations  in  order  to  achieve  a  final  permutation  for  which  Z  is  close  to  its  global 
minimum.  The  need  for  multiple  starts  increases  the  computational  effort  many-fold. 

1 .  First  Example 

A  16x16  example  from  Ref.  33  was  solved  with  the  gradient  algorithm  for  moment 
compression.  In  this  example,  the  16  most  frequently  occurring  non-trivial  words  have  been 
extracted  from  a  long  conversation.  The  input  matrix,  A-l ,  is  shown  in  Fig.  27.  A  “1  ”  is  placed 
in  row  i  and  column  j  if  words  i  and  j  have  coincidentally  occurred  in  two  or  more  sentences, 
and  a  “0"  is  placed  there  otherwise. 
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FIGURE  27.  Initial  Word  Relationship  Array,  Matrix  A-l 


Since  the  input  matrix  is  symmetric,  the  problems  of  choosing  permutations  to 
minimize  the  row  or  column  moments  are  identical.  It  sufficed  to  find  the  optimal  column 
permutation,  and  to  use  this  permutation  on  both  the  rows  and  columns  of  the  matrix.  The 
problem  of  row  moment  minimization  was  solved  40  times,  each  time  starting  from  a 
randomly  chosen  permutation  of  the  columns.  Two  solutions  are  taken  as  identical  if  they 
differ  merely  bv  reversal  of  the  order  of  the  16  words. 

The  results  were  as  follows.  Nine  of  the  40  starting-points  led  to  the  final  matrix,  A-2, 
shown  in  Fig.  28,  with  the  lowest  ME.  An  additional  1 1  of  the  40  starting-points  led  to  a  final 
matrix  (with  the  same  ME)  which  differed  from  A-2  only  by  interchanging  of  the  variables 
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FIGURE  28.  Reordered  Word  Relationship  Arrays 
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bed  and  “lips”.  Since  these  two  variables  have  identical  rows  and  columns,  it  is  under¬ 
standable  why  ambiguity  arises  about  their  ordering.9 

Inspection  of  Fig.  28  shows  that  the  algorithm  has  partitioned  the  16  words  into  three 
subjects:  vacation,  sex,  and  sports.  The  transitional  words  between  these  three  topics  of 
conversation  are  evidently  hotel,  hot,  court,  and  leg. 

All  remaining  20  starting-points  led  to  a  ME  which  was  at  least  1 1  percent  oighev.10 
Inspection  of  the  next-to-best  permutation  (the  one  with  an  ME  1 1  percent  higher  than  that 
of  matrix  A-2)  showed  that  the  gradient  algorithm  converged  to  the  wrong  local  minimum  of 
Z(X),  in  which  only  one  of  the  three  topics  of  conversation  (sports)  is  clearly  identified. 

It  is  believed  that  Matrix  A-2  (and  the  variations  obtained  by  permuting  identical  rows) 
achieves  the  global  optimum  of  Z(X),  although  this  is  not  certain.  It  must  be  recalled,  however, 
that  the  primary  goal  is  the  discovery  of  informative  patterns,  not  rigorous  optimization  of  the 
ME.  For  example,  the  rearrangement  proposed  by  Giuliano,  Matrix  A-3  in  Fig.  28,  is  just  as 
pleasing  as  A-2,  even  though  its  ME  is  not  as  good.  The  point  here  is  two-fold:  (1)  data 
rearrangements  with  near-optimal  ME  may  be  as  pleasing  as  those  with  optimal  MEs;  (2) 
ME-optimizarion  algorithms  can  fail  to  identify  all  informative  patterns,  especially  patterns 
which  are  not  local  optima  for  the  ME. 1 1 

'Hie  Moment  Compression  Algorithm  is  considered  to  have  worked  properly  on  this 
example  because  it  produced  a  pleasing  pattern.  The  sensitivity  of  the  gradient  algorithm  to 
the  starting  point  was  not  a  serious  problem,  for  20  of  the  40  starts  led  to  a  good  answer. 
Note,  however,  that  a  mere  1 1  percent  degradation  in  the  ME  led  to  a  seriously  degraded 
pattern. 

The  main  criticism  of  the  gradient  algorithm  for  Moment  Compression  is  its  excessive 
computation  time.  Each  usage  of  the  algorithm  required  3  to  7  gradient  steps  (i.e.,  solutions  of 
3  to  7  linear  assignment  problems)  at  about  2  seconds  per  step.  The  algorithm  therefore 
required  about  10  seconds  per  starting  point.12  Since  40  starting  points  were  chosen  at 
random,  to  ensure  high  confidence  in  achieving  a  global  rather  than  local  optimum,1 3  7 
minutes 14  were  required  to  solve  this  problem. 

9.  It  should  bo  pointed  out  that  this  example  exhibits  considerable  degeneracy.  Examination  of  Eig.  33  reveals  that 
rows  7-9,  rows  10-11,  and  rows  14-16  are  identical.  The  ME  will  be  invariant  under  permutations  of  identical  rows.  In 
addition,  rows  t-4  and  rows  5-6  are  nearly  identical;  the  ME'  will  undergo  only  minor  changes  if  these  are  interchanged. 

10,  The  ME  Is  here  defined  as  the  root-mean-square  row  moment. 


1 1 .  Since  the  sterling  points  never  led  to  A-3,  A-3  probably  is .  >st  a  local  minimum  for  Zi  X  i. 

13.  By  contrast,  the  Bond  Energy  Algorithm  requires  only  a  few  seconds  pet  slatting  point. 

13.  By  contrast,  the  Bond  Energy  Algorithm  is  rather  m sensitive  to  the  slatting  point,  so  that  fewer  starting  points  (at 
most  16,  and  ptobably  much  less)  need  be  explored. 

14.  On  CDO  1604. 
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It  may  be  possible  to  cut  the  running  time  by  a  factor  of  2  or  4,  since  a  “good”  starting 
point  may  require  fewer  gradient  steps  than  a  randomized  one,  and  by  using  much  fewer  than 
40  starts.  Nevertheless,  the  running  time  for  a  16x16  problem  (on  the  order  of  minutes,  and 
doubled  if  the  matrx  is  not  symmetric)  is  disappointing  when  contrasted  with  the  running 
time  for  the  Bond  Energy  Algorithm.  Consequently  the  gradient  algorithm  for  Moment 
Compression  Algorithm  is  probably  impractical  for  problems  larger  than  25x25  or  so.  It  works 
well,  but  is  too  slow. 

2.  Second  Example1  s 

A  second  example  was  run  in  order  to  test  the  ability  of  the  gradient  algorithm  to 
generate  clumps  of  large  numbers  when  the  matrix  elements  were  not  restricted  to  0  or  1 .  The 
initial  matrix  is  the  10x10  matrix  denoted  as  B-l  in  Fig.  29.  Since  B-l  is  symmetric,  it  sufficed 
to  find  the  optimal  row  permutation,  and  to  use  this  permutation  on  both  the  rows  and 
columns  of  B-l. 
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FIGURE  29.  Initial  Similarity  Matrix  B-l 


The  gradient  algorithm  was  used  with  60  randomly  chosen  starting  points.  Four 
distinct  MEs14  were  obtained,  with  values  1.846,  1.987,  1.988,  and  2.314.  The  frequency  of 


13.  This  is  the  «ame  example  as  in  Figs.  7  and  8, 

16.  The  ME  is  here  defined  as  the  tool-mean  square  tow  moment -arm, 
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these  MEs  were,  respectively,  36,  12,  8,  and  4  (sum  of  60).  The  matrices  corresponding  to  the 
four  MEs  (the  best  four  permutations17)  were  respectively,  B-2,  B-3,  B-4,  and  B-5  and  are 
shown  in  Fig.  30.  These  four  MEs  are  fairly  close  to  one  another  and  it  is  apparent  that  the 
patterns  in  B-2,  B-3,  B-4,  and  B-5  are  essentially  equivalent  and  equally  pleasing;  all  four 
matrices  succeed  in  identifying  a  2x2  block  of  large  numbers  (variables  H  and  D),  a  4x4  block 
(variables  A,  B,  E,  and  1)  and  a  4x4  block  (variables  C,  F,  G,  and  J).  The  exact  order  of  the 
blocks,  and  of  variables  within  each  block  differ,  but  one  would  be  indifferent  to  such 
unimportant  differences  (i.e.,  to  the  arrangement  of  the  stray  l’s)  since  the  identification  of 
the  blocks  of  primary  blocks  is  what  is  significant. 

The  computation  time  for  this  problem  was  about  3  seconds  per  starting  point.  Tire  60 
starting  points  consumed  about  3  minutes  total  computation  time1 8. 

Summarizing,  the  gradient  algorithm  applied  to  this  problem  succeeded,  in  all  60  of  the 
starts,  in  identifying  the  major  variable  blocks  and  produced  informative  patterns.  The  four 
best  permutations  produced  nearly  equal  MEs  and  equally  informative  patterns,  without  any 
obvious  way  of  choosing  among  them.  The  algorithm  is  considered  successful,  but  rather  slow, 
for  this  problem.  It  correctly  “factored”  the  main  variables,  but  was  very  time-consuming 
compared  with  the  Bond  Energy  Algorithm. 


17.  Two  permutation!  were  eomtileieti  equivalent  tl  each  coulj  be  obtained  Horn  the  other  by  merely  revervtrtg  the 

outer  of  the  10  row*  ami  columns. 

18.  Go  CPC  1604  computer. 
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FIGURE  30.  Reordered  Similarity  Matrices 
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FORMULATION  OF  THE  BOND  ENERGY  ME  OPTIMIZATION 
AS  TWO  QUADRATIC  ASSIGNMENT  PROBLEMS 


The  purpose  of  this  Appendix  is  to  show  how  the  problem  of  Mi:  maximization  can  be 
rigorously  formulated  and  solved  as  two  Quadratic  Assignment  Problems  (QAPs).  This  form¬ 
ulation  is  presented  only  for  theoretical  interest,  because  published  algorithms  (Refs.  21,  22, 
23)  which  find  truly  optimal  solutions  to  QAPs  are  too  time  consuming  to  be  practical  for 
large  problems.1  While  approximate  algorithms  have  been  published'  (Refs.  24,  25,  26,  27) 
which  find  near-optimal  solutions  to  QAPs,  it  was  not  believed  worthwhile  to  explore  any  of 
them,  because  none  exploited  the  nearest-neighbor  feature  of  the  function  being  optimized. 
Only  the  sequential-selection  approximate  algorithm  described  in  this  paper  exploits  the 
nearest-neighbor  feature,  and  this  latter  algorithm  is  believed  to  be  much  faster,  more 
convenient,  and  just  as  satisfactory3  as  the  published  approximate  QAP  algorithms. 

Suppose  the  original  non-negative  matrix  t  ay  1  is  M  x  N  with  horizontal  and  vertical 
bond  energies  contributing  to  the  ME.  The  ME  then  consists  of  the  sum  of  two  terms,  namely, 
the  row  bond  energies  plus  the  column  bond  energies.  Two  optimization  problems  must  be 
solved  for  ME  maximization.  One  seeks  a  permutation  of  the  columns  of  |  ay  1  which  maximizes 
the  row  bond  energy,  the  other  seeks  a  permutation  of  the  rows  of  [ay !  which  maximizes  the 
column  bond  energy.  These  two  optimization  problems  can  be  carried  out  independently  of 
each  other.  When  both  are  completed,  the  optimal  permutations  of  both  rows  and  columns  are 
known. 


The  two  optimization  problems  are  mathematically  equivalent.  Only  the  problem  of 
maximizing  the  row  bond  energy  is  presented  here.  This  problem  requires  selection  of  a 
permutation  ir  =  [*(1),  jt(2),  ....  tr(M)|  of  the  integers  [1,2, ....  M|  which  maximizes 


N  ( 


bil  bi2 


M 


j  =  2 


bU 


Ou- 


+  b: 


1] 


+  b 


iM  °i,  M  -  I 


( A-I ) 


1.  Computer  Umes  on  the  order  of  one  or  seve>al  minutes  are  required  for  15x13  matrices.  uml  use  as  the  fourth  ami 
fifth  power  of  the  matrix  sine. 

2.  An  extensive  bibliography  ts  contained  in  Ref.  26. 

J.  The  satisfaction  with  the  Bond  Energy  Algorithm  is  not  based  primarily  on  how  close  it  comes  lo  achieving  me 
global  optimum  in  Eq.  (A-2)  but  rather  ou  the  pleasing  patterns  of  clumps  which  n  produces. 


Preceding  page  blank 
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The  term  within  braces  is  twice  the  bond  energy  for  the  i**1  row  of  [by],  where  [by]  =  [ay^] 
denotes  the  matrix  whose  j*k  column  is  the  *0')^  column  of  [ay]  .  The  mathematical 
problem  may  be  rewritten  as 

M  M 

max  £  £  Qjkir(j)jr(k)  (A‘2) 

n  J- 1  k  =  1 


Qjkrs  ,£j  air  ais  L^k,  j- 1  +  ^k,j  +  l] 

1  <j,  k,  r,  s  <M 

The  maximization  in  Eq.  (A-2)  is  taken  over  all  M!  possible  permutations.  This  type  of 
maximization  is  known  as  a  quadratic  assignment  problem  because  of  the  double  appearance 
of  ir  in  the  maximand.  As  previously  noted,  published  algorithms  exist  for  finding  both 
optimal  and  near-optimal  solutions  to  Eq.  (A-2). 

INTERPRETATION  OF  ME  OPTIMIZATION  AS  TWO  TRAVELING  SALESMAN  PROBLEMS 


The  quadratic  assignment  problem  formula' cd  in  the  previous  section  is  actually  a 
special  type  called  the  open-loop  traveling  salesman  problem.  Let 

N 

^rs  =  £  airais  "  ^sr 
i  =  1 

denote  the  scalar  product  of  the  r^  and  s”1  columns  of  [ay].  Then,  the  maximization  in  (A-2) 
is  equivalent  to 

M-  1 

max  £  ^ir(j)ir(j  +  1) .  (A-4) 

v  J  “  * 

If  one  interprets  drs  as  the  distance  4  from  city  r  to  city  s,  the  problem  in  Eq.  (A-4)  is  to  find 
the  salesman's  tour  [from  city  ir(l)  to  w(2)  ..  to  city  »(M)j  of  the  M  cities  whicn  has  the 
longest  distance.5  Note  that  the  tour  origin  is  arbitrary  and  that  the  salesman  is  not  required  to 
return  to  his  origin.  This  tour  is  therefore  called  open  loop. 

4.  It  necessary,  a  large  positive  constant  can  be  added  to  nil  dt  in  order  to  nuke  them  positive. 

5.  Subtraction  of  every  d[s  from  a  large  positive  constant  kadi  to  an  equivalent  problem  of  minimizing  the  tour 

length. 
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PROOF  THAT  THE  BOND  ENERGY  SUBOPTIMAL  ALGORITHM  WILL  PRODUCE 
BLOCK  FACTORED  FORM  IF  IT  IS  POSSIBLE  TO  DO  SO 
BY  ROW  AND  COLUMN  PERMUTATIONS 


The  purpose  of  this  appendix  is  to  prove  that  the  sequential  selection  bond  energy 
algorithm  will  put  a  matrix  into  block  factored  form  if  it  is  possible  to  do  so  by  row  and 
column  permutations. 

DEFINITION  1 

A  non-negative  matrix  A  whose  elements  ay  relate  row  entity  i  to  column  entity  j  is 
called  block  factorable  if  the  row  entities  can  be  decomposed  into  q  disjoint  subsets 
R  |  ,Ri.  •  •  •  Rq>  and  the  column  entities  decomposed  into  q  disjoint  subsets  C j  ,C->. . . .  Cq  with 
the  properties: 

(!)  If  entity  i  e  Rfl.  then  ay  =  0 

unless  entity  j  e  Ca,  l<a<q 
and  if  entity  j  e  C'rt,  then  ay  =  0 

unless  entity  i  e  Rfl,  1  <  a  <  q 

(2)  For  each  a.  the  submatrix  j  [ay  ]  .  ieRfl,  jet’a^ 

cannot  be  further  decomposed. 

That  is,  A  can  be  factored  into  q  blocks  if  the  row  entities  and  column  entities  can  each  be 
partitioned  into  q  subsets  such  that:  ( l )  entities  in  one  row  subset  interact  only  with  entities  in 
the  corresponding  column  subset  and  ( 2)  it  is  impossible  to  decompose  the  subsets  further. 

DEFINITION  2 

A  block  factorable  matrix  is  said  to  be  in  block  factored  form1  when  all  the  row 
entities  contained  in  each  Ra  lie  together  on  the  vertical  axis  of  the  matrix  and  all  the  column 


t.  Figure  4  show*  a  matrix  in  block  factored  form. 

Preceding  page  blank 
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entities  contained  in  each  Cfl  lie  together  on  the  horizontal  axis  of  the  matrix.  Clearly,  the 
matrix  A  is  block  factorable  if,  and  only  if,  it  can  be  put  into  block  factored  form  by  row  and 
column  permutations. 

LEMMA  1 

Assume  A  is  block  factorable.  If  row  entity  i  of  matrix  A  is  contained  in  RQ  and  row 
entity  j  is  contained  in  with  a^, then  the  scalar  product  of  row  i  with  row  j  vanishes. 

Proof 

For  any  entity  k,  2^=0  unless  keCa  and  ajjc=0  unless  keC^.  Therefore, 
aik  f°r  k  s‘nce  a  anc*  Cp  are  disjoint. 

LEMMA  2 

Assume  A  is  block  factorable  and  select  any  Rfl  which  contains  two  or  more  rows.  No 
matter  how  Rfl  is  split  into  two  distinct  subsets,  it  is  always  possible  to  choose  one  row 
from  each  subset  such  that  the  scalar  product  of  the  two  rows  is  positive. 

Proof 

If  such  a  choice  cannot  be  made,  then  the  submatrix  |(ay],  ieRfl,  jeCa  }  is 
decomposible,  violating  Definition  1. 

THEOREM 

If  A  is  block  factorable,  then  the  sequential  selection  algorithm  will  put  the  matrix  into 
block  factored  form,  and  will  do  so  by  building  one  block  at  a  time. 


If  the  first  row  laid  down  came  from  (say)  Rj,  then  the  next  row  to  be  laid 
down  will  be  one  of  the  remaining  M-l  rows  with  the  greatest  scalar  product 
with  the  first.  Since  (by  Lemma  I)  all  the  rows  not  contained  in  Rj  have 
vanishing  scalar  products  with  the  first  row,  and  since  at  least  one  (by  Lemma 
2)  of  the  as  yet  unplaced  rows  from  Rj  (if  any  others  exist)  has  a  positive 
scalar  product,  then  the  second  row  to  be  laid  down  will  come  from  Rj.  By 


86 


repeating  this  reasoning  it  is  clear  that  all  the  rows  from  Rj  are  laid  down 
before  any  other  rows  are  laid  down.  More  generally,  one  subset,  Ra,  of  row 
entities  at  a  time  is  laid  down  and  all  the  rows  contained  in  each  Rfl  lie  together 
in  the  matrix. 

Identical  reasoning  can  be  applied  to  show  that  the  columns  are  also  laid  down 
with  all  the  columns  in  each  Ca  lying  together.  Therefore,  by  Definition  2  the 
matrix  will  be  put  in  block  factored  form. 
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APPENDIX  C 


PROOF  THAT  THE  MOMENT  COMPRESSION  ALGORITHM  WILL  PRODUCE 
BLOCK  FACTORED  FORM  IF  IT  IS  POSSIBLE  TO  DO  SO 
BY  ROW  AND  COLUMN  PERMUTATIONS 
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PROOF  THAT  THE  MOMENT  COMPRESSION  ALGORITHM  WILL  PRODUCE 
BLOCK  FACTORED  FORM  IF  IT  IS  POSSIBLE  TO  DO  SO 
BY  ROW  AND  COLUMN  PERMUTATIONS 


The  assertion  here  is  that  the  minimum  of  the  Mr.1  of  a  matrix  which  can  be  placed  in 
block  form  via  row  and  column  permutations  occurs  when  the  matrix  is  in  block  form,  and 
does  not  occur  when  rows  (or  columns)  of  one  block  are  separated  by  rows  (or  columns)  of 
another  block.  Consequently,  rigorous  minimization  of  the  ME  must  put  the  matrix  into  block 
form.  Since  the  gradient  algorithm  for  moment  compression  will  find  a  global  rather  than  local 
minimum  of  the  ME,  if  sufficiently  many  starting  points  are  used,  it  follows  that  the  gradient 
algorithm  will  put  a  matrix  into  block  form  if  this  is  possible. 

It  suffices  to  examine  how  column  permutations  can  minimize  Zj  rj.  The  basic  idea  of 
the  proof  is  that  if  the  columns  from  one  block  are  separated  by  columns  from  another  block, 
then  removal  of  the  extraneous  columns,  reuniting  the  columns  from  the  first  block,  and 
reinsertion  (at  one  side)  of  the  removed  columns  will  strictly  reduce  Zj  q,  hence  reduce  the 
ME.  Thus,  the  ME  is  at  its  minimum  only  if  columns  from  the  same  block  are  contiguous. 

An  example  is  provided  by  Fig.  C-I  which  shows  the  5  left-most  columns  of  a  matrix. 
At  least  one  X  in  each  column  is  positive.  Columns  A,  C,  E  form  a  block;  no  column  to  the 
right  of  E  lies  in  this  block;  and  columns  B  and  D  are  from  other  blocks. 

The  following  theorem  shows  that  if  column  D  is  moved  out  to  the  right  of  the  block 
(producing  the  column  order  A,B,C,E,D),  then  ZjTj  will  decrease.  Similarly,  movement  of 
column  B  to  the  right  of  the  block  (producing  the  column  order  A,C,E,B,D)  will  reduce  Zjq 
further. 


1. 

N  M 

ME  =  £  rj  +  L 
i=  l  )=  I 


moment  arm  for  row  i 
moment  arm  for  column  j 


Preceding  page  blank 
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COLUMN;  A  B  C  D  E 

X  0  X  0  X 

X  0  X  0  X 

X  0  X  0  X 

0X0X0 
0X0X0 
0X0X0 
0X0X0 

FIGURE  C-1.  Sample  Matrix 


The  general  procedure  is  to  identify  the  left-most  column  block  whose  columns  are  not 
placed  contiguously,2  and  to  move  the  right-most  extraneous  column3  from  the  midst  of  the 
block  to  the  immediate  right  of  the  block.  Repetition  of  this  procedure  produces  a  column 
ordering  which  places  columns  from  the  same  block  in  contiguous  positions.  Since  the 
procedure  leads  to  strict  decreases  in  Xjrj,  it  shows  that  the  ME  achieves  its  global  ME  only 
when  the  matrix  is  put  into  block  factored  form. 

The  theorem  which  follows  shows  that  each  rj  is  decreased  if  the  zeros  which  are 
interior  to  a  row  are  moved  to  an  edge  of  the  block,  and  is  unchanged  if  a  zero  at  one  block 
edge  is  moved  to  the  other  block  edge.  Thus  moving  columns  B  and  D  to  the  right  of  E  will 
reduce 

3 


because  the  first  3  rows  have  an  interior  (or  possibly  left  edge)  zero  at  column  B,  and  an 
interior  (or  possibly  right  edge)  zero  at  column  D.  Similarly, 


is  also  decreased  by  such  a  transfer  because  columns  C  and  E  provide  interior  zeros  to  rows 
4-7. 


2.  Initially,  this  is  block  A.C.E. 

3.  Initially,  this  is  column  O. 


THEOREM 


N  N 

Let  Wj>0,  £  W:  =  1 ,  j  =  £  jWj 

j=l  j=l  J 

N 

S  =  £  wi  C  J  ~  i  ]  ^  ^  moment  for  the  vector  W. 
j  =  1  J 

Suppose  for  some  k,  2<k<N- 1 ,  W^O.  Let  *  refer  to  a  rearrangement  whereby  Wk  has 
been  moved  to  the  extreme  right,  and  the  vector  then  closed  up: 


* 

_  I 

Wj 

1  <j<k-  1 

wi 

Wj  +  1 

k<j<N  -  1 

(C-l) 

0 

j  =  N 

J 

III 

• 

II 

- * 

=  j  -  j3  where  j3  = 

N 

E  ^ 

J  =  k  +  1  J 

(C-2) 

a 

N 

s 

=  I 

j  = 

'  w* 

J 

r  .  — -i  £ 

1_J-J  J  =  the  moment  for  the  vector  W*. 

(C*3) 

Then  S*  <  S  with  equality  if  and  only  if  WR  is  an  edge  zero  (that  is,  if  W~0  for  all  j<k.  or  if 
Wj=0  for  all  j>k). 


Proof: 


Set  j  =  E  +  F  where 


k-  1  k-  1 

E3  £  W:j<(k-1)  £  w:  3  (1 -0)(k -  1) 

i=i  ’  i-i 


N  N 

F  =  £  ws  j  Xk  +  I)  E  wi  = 

j-krl  J  j3k+l  ' 


(C-4) 


(C-5) 
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Insertion  of  Eqs.  (C-l ,  C-2)  into  Eq.  (C-3)  obtains 

*  N 

S  -S-0(l-0)  +  2  V  W;(j-j) 
j=T+ 1  J 

=  0(1-0)  +  2  0j  -  2F  =  0(1-0)  +  20E  +  2(0-  I)F 

Insertion  of  Eqs.  (C-4,  C-5)  produces,  because  0  <  0  <  1 ,  the  result 
*  * 

S  -  S  < -  30  ( 1  -  /3)  <  0  with  S  =  S  only  when  0  is  0  or  1,  which  occurs  only  if  w^  is  an 
exterior  zero.  Thus  S*  =  S  only  if  is  an  exterior  zero.  The  converse  is  easily  proved.  QED. 

Note  that  with  the  choice 


/  N 

Wj *  au/  t  *im 

/  m  =  1 

wefindS  =  rj  =  moment  for  itfl  row  of  {ay}.  The  theorem  therefore  states  that  the  *  - 
rearrangement  (namely  removal  of  an  interior  zero  from  the  i^  row  of  [ay])  will  strictly 
reduce  q  unless  the  zero  lies  at  an  edge. 


94 


APPENDIX  D 


THE  BOND  ENERGY  COMPUTER  PROGRAM 
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THE  BOND  ENERGY  COMPUTER  PROGRAM 


A.  OPERATION  OF  THE  PROGRAM 

The  computer  program  for  the  Bond  Energy  Algorithm  consists  of  two  separate  parts. 
The  first  part  of  the  program  reorders  the  columns  while  ihe  second  part  reorders  the  rows. 
Figure  D-l  shows  the  essential  program  logic  for  selection  and  laydown  of  the  rows  to  obtain  a 
new  order  with  a  large  NME.  The  logic  for  the  column  selection  is  identical.  It  was  found  that 
in  order  to  examine  a  number  of  local  minima  it  was  necessary  to  initiate  tnc  program  at 
several  starting  rows  (or  columns).  However,  as  pointed  out  in  Chapter  1  of  Part  11.  almost  all 
starting  points  (rows)  resulted  in  a  “good”  solution. 

B.  CARD  INPUT  FORMAT 

The  input  format  that  :.s  described  here  is  for  arrays  with  integer  elements.  The  only 
change  that  would  have  to  b  :  made  to  accommodate  decimal  entries  is  in  the  input  and  output 
formats  for  the  initial  and  reordered  arrays. 

nat (415) 

number  of  rows  in  the  matrix 
number  of  columns  in  the  matrix 
an  increment  to  determine  the  starting  columns  or  rows. 

If  IFZI  =  5  then  the  algorithm  is  run  K  times  beginning 
with  row  1  then  row  6,  and  continuing  in  increments  of 
5  until  K*IFZl  +1  >  NN  or  MM 
0  or  blank  if  the  input  array  is  not  symmetric 
1  if  input  array  is  symmetric 

b.  Cards  2  through  MM  +  1  Format  (8011) 

(NA(I,J),  J  =  1,  NN)  is  the  l^1  row  of  the  mput  matrix.  This  card 
is  repeated  for  all  MM  rows. 


a. 


Card  1  Forr 
MM 
NN 
1FZ1 


IFSYM  = 
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FIGURE  D-l.  Flow  Chart  for  Sequential  Row  (or  Column)  Selection  ond  Laydown 

for  Bond  Energy  Algorithm 
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C.  PROGRAM  OUTPUT 

The  computer  program  GROUP  2  consists  of  the  following  output  information: 

( 1 )  A  printout  of  the  input  cards. 

(2)  A  printout  of  the  new  row  and  column  orderings  at  each  step  in  the 
sequential  laydown  procedure. 

(3)  A  printout  of  the  final  reordered  matrix. 

(4)  A  printout  of  the  final  horizontal  and  vertical  MEs. 


f 
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BOND  ENERGY  ALGORITHM  PROCRAM 


COMMON  MM.NN 

COMMON  NA(90,90>,k<9°»  »L<9°)*ME<«0)  , NPOSl 9O ) » NS* 

set  i npe x  a  •  *.  h  ■  1 

IDIM  ■  90 

Bf*D  lOO.HM.NN,  JF21,I|fSVM 
irtirsvM.EO,*0)  irsvM.p 
P*INTlno.MM4NN,IE2l.IFSYM 
100  F0Rm*TMI5» 

CO  1  I ■ 1 #  MM 

HI*D  101*(NA(I4IJ)4Jb14NN} 

DO  12  J.1.NN 

IriNilliJi .£Q.-0 )  N*(I,JJ«0 

12  continue 

P9INT  102* (N*(I,JJ»J«1,NNJ 
1  CONTINUE 
1 91  F OBm* Y  t  gftll > 

102  E OflM/iT <  1 X,  J9I2  ) 


DO  22  KKn.NN.JEIl 

DO  2  In. mm 

2  *H).i 

DO  3  JkI.NN 

3  „(J>  •  TOJM  •  (J«l) 

K  1 1 ) «KK 
K(KK)*i 
ITEMPl  .  L<1) 

111)  ■  L I  KK ) 

UKK)  ■  ITEMPl 

NTICK«0  t  MTICK«0 
M|1»0«mf?«0 

Nl»NN»i 


03  150  NCOUNTn.Nl 
N(1  bNCOijNT  *1 

DO  200  JkNCI.NN 
M|l J)»0 
N90S I J ) • 0 

DO  310  NIO.NCOUNT 

nuM.iooo 

DO  211  j«l.MM 

G  •  K(Xt  *  LU) 

M  •  K«I)  •  L<N) 

M  •  K(t,  *  U*Ui> 

lEI N ,CO . 0 ISO  TO  212 
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ST’iyv" 


If  INaiGI.EQ.0)  so  to  213 
If  I NA I M) , E  0. 0  )  SO  TO  ?12 

NJUH  ■  nSUH  *  NAIM)  *  N  a  <  a  > 

212  If (N.EO.NcOUNTiSO  to  211 
If  INAIM). E0.°>  SO  TO  211 

NJUM  ■  SSU*»  *  NAIM)  •  N*(S) 

213  IfIN.EQ.fl.Ofl.N.EQ.NcOUNT)  SO  To  211 
If  INAIM). EQ.O)  GO  TO  211 

N5UM  a  nSU“  »  NAIM)  *  NAIM) 

211  CONTINUE 

IflNSUM.LT.MEI J) ISO  TO  210 
IflNCOUNT.fO.il  GO  TO  214 
IflNSUM.NE.MEIU)>  GO  TO  214 
If IN.EC.O.oa.N.EO.NcOUNT)  SO  TO  214 
GO  TO  ?iA 

214  CONTINUE 

ME | J1BNSUM 

Nf OS t J ) bN 
210  CONTINUE 

200  CONTINUE 


M£MAXBMf( NCI > 

N»maxbNPOSINcD 
IMA  XiNCi 
NC2»NCOl|NT*2 
00  220  I«NC2,NN 
NSHsl 

If  IMEID.LT.m|maX)GO  TO  220 

IflMEID.E8.MEMA*>  OKI  £82l< I*  1»A* ) 

If  I  NSW . NE . 1 )  GO  TO  220 

IMAX.J 

M£M a  XbME I  X ) 

NPMAXbNPOSII) 

220  CONTINUE 

HElBMCi *M|M*X-1000 
If|NPMtx,NE.NC0UNT)QO  TO  l4S 
NNl«L< IMA*) 

UlMA*)iLINPMAX»l) 

LI NPMAX*1 I *NNl 
GO  TO  1J0 

140  If INPMax.NE.O)  NTIC*"NTXC**1 
NPIbNPmaXbI 
NSAVEbuIMAX) 

NP2»Nfi«l 

FOR  145  I«IMA*,NPa»l 

L 1 1 )  *L  ( x»x  ) 

149  CONTINUE 

LINPDbNSAVE 

190  *«INT  109* l L I  X  )/XD!M*l«  X«1*NN) 


lfllfSVM.CO.O)  GO  TO  191 

DO  1 52  1*1 »NN 
192  MX)*III)  /  XDIM  *  1 
M(2*HCi 
MfXCN»NTICN 
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'  "  ' •  •*«  VR^'rap  -/;  r- 1 4  3  7 


*  v-  'a  ^ywv.v».T  ?,Tf»^  y?,',,ivr-~-Z'Vr:-fi'\’*^Tfi 


GO  TO  S$1 


151  continue 

105  FOBMATUOIS) 

Hl»MH-5 

DO  550  MC0UNT«1,M1 
MC1«HCOuNT*l 

DO  400  I«MC1,MM 
f I > «° 

NPOS  ( I )  jT* 

DO  410  m*0#mcOUNT 
MSUMslOOO 

DO  411  J * 1 » N N 

G  »  K  < I )  «  U(d) 

H  »  K<K)  »  L ( J ) 

N  »  K<“*1 I  ♦  U  <  J  > 

lPfH.EO.")GO  TO  412 
IP  (N*(G)  .EQ.O)  go  TO  413 
IP  (NAfMl .EQ.O)  60  TO  412 
M$UM  ■  mSUm  ♦  NA(H)  *  N*(G) 

412  Ir(M,EO.*cOUNT>GO  TO  413 
IP  (Na<N> .CQ.Ol  60  TO  411 
MJUM  ■  m5Um  ♦  NAfNl  *  NA(6) 

413  IP(M,EQ. 0 .PR.M.CO.MCOUNT)  GO  TO  411 
IF  (NAfHl.EQ.Ol  60  TO  4ll 

M$UH  ■  mSUX  -  N  A  f  H  1  *  MAIN) 

4ll  CONTINUE 

iFlHSUM.lT.xemiQO  TO  410 
lPfMcOUNT.EQ.il  60  TO  414 
lF(MSUM,Nf .Mffl)  1  GO  TO  414 
iPfM.CO.n .OR.M.EQ.MCOUNTl  G»  TO  41« 
60  TO  «ln 

414  CONTINUE 
MC«I)»MSUM 
NPOS  ( X 1  ■** 

4l 0  CONTINUE 

400  CONTINUE 


M|MAX»Mf  t MCI • 

NPMAXbNPOSI MCll 
JMAXaMCl 
Mca«MC0uNT*2 
DO  4^0  J-MC2.NM 

lPfME(Il.lT.M|MAX)CO  TO  420 
IPIMEfli  .lO.MEMAlf)  C*U  E*tJ<X»d**X> 

4  * 


XPtNSW.NC 
J*AX«X 
M£MAX*MEf  Jl 


)  00  TO  420 
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~  1 1  K~ 


NPNAX«NPOS<I> 

420  CONTINUE 

mE2»HE2*m|max-1000 
If(NPMAK.NE.HCOUNT>GO  TO  J40 
MMl«K(JMAX J 

R(  JMAX)»*<NPMAX*i  ) 
K(NPMAX*il»MRi 

GO  to  S5n 

340  XPI NPMax  ,N| ,  0 >  MTXCK»MTIC**i 
MR1«NP>"AX*1 


345 

350 

351 

352 


5OO 

22 

502 

501 

SO) 


MSAVEM  f  JMAX) 

MP2«*Pi +1 

FOP  345  ,j*jMAX„MP2,l 
«IJ)  > 

CONTINUE 

k|mPH«m$aVE 

PRINT  1 0 5 #  <  K t  J » . J*i «  mM} 

DO  352  5  •  1,NN 

l(G>  •  L<a>/IDT«  •  1 

PRINT  50s,  <U J>» Jil.NM 

DO  500  x«l,NM  « 

PRINT  501  ,K<I>  ,<NA<K<I>  jUU> 

PRINT^S 03  «NTlCX*NTICK«*<El|Ni2 

F0Bhat'(//2I5»5X**«1***15‘5|‘**mC2***1,//> 
FOPNATf ix,l2,  ?X»|7I</> 

FORMAT ( ini ,//5X, 3712// > 

END 


subroutine  eq2HH*xm> 

COMMON  MM.NN 

COMMON  NAlfO ,»0)»R19°)»L 

set  index  a  •  x»  *  •  t 


(t0),MC<90>,NB0S(*0)*N5N 


NSN«0 

XflNPOS<IX« ,E8.NROS(IR>»aO  TO  701 
RETURN 


701  DO  702  I«liHM 

a  ■  mu  •  i<n> 

N  ■  K  1 1 |  •  l(XM> 

IP  (NA«Q>  .N£.NA<MJ) 

702  CONTINUE 
return 


HO  to  7«» 


70)  N|UM1»0 

N|UN2*0 
DO  704  X»l.«* 
a  ■  mix)  •  uii> 

M  •  KID  •  UlXM> 

NIUMl  ■  NfUMi  ♦  NA»G» 
704  N*U»2  •  NfUMJ  ♦  NM*41 

XPIN*UMi.LT.NlUM2)NSM«l 

return 

END 


_ 

% 


$ 


I 

4# 


I 

*& 


© 


l 


(■ 

mu 


l 

9t 


‘i 

*& 


i 

l 

i 

l 

l 

l 


SUBROUTINE  C02J(JJ»JM) 

COMMON  mm.nn 

COMMON  NA(«0.9<>)  <K<9C),ittO)#MC<eO  ),NP0S<  «0)#NSN 

set  inpe*  o  •  *,  h  •  l 
NS*»0 

ir<NPOS( jjj ,eo.npos< jm>go  to  701 

RETURN 

701  Do  702  1*1, NN 

3  ■  K « JJ  >  ♦  L  ( I  > 

M  «  MJM>  •  L  <  T  > 

1 1  <  N  A  (  Q  )  .NE.NMMJ)  GO  TO  70} 

702  CONTINUE 

return 


7qJ  N5UMl*(t 

N|UM2*0 
DO  704  I* t » NN 
G  *  K  (  j  j )  «  LIU 
"  ■  KjjMJ  *  l  <  1 ) 

NSUMl  ■  KSUmI  *  N*(5) 
7O4  NSUM2  ■  w S U m 2  •  NA(H) 

If <NSU"i.LT.NSUM2)NS“*t 
return 

END 


subroutine  eoi 
c  temporary  sub 
print  «oe 

»09  f ORm*T ( «F NT£R  £01* J 
End 
End 
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THE  MOMENT  ORDERING  COMPUTER  PROGRAM 

A.  INTRODUCTION 

Two  versions  of  the  program  are  available.  One,  which  handles  arrays  of  size  up  to 
100x16,  reads  the  arrays  directly  from  cards.  The  second,  which  is  identical  to  the  first  except 
in  its  input-output  procedures,  can  handle  arraysupto  100x100 and  reads  its  arrays  from  tape.1 
Because  the  two  are  so  similar,  only  the  first  will  be  presented  here;  most  of  the  discussion, 
however,  applies  equally  to  both  versions.  A  flow  chart  showing  the  main  program  logic  is 
given  in  Fig.  E-l.  The  following  section  provides  instructions  for  the  use  of  the  program; 
Section  C  describes  the  program  output. 

B.  INSTRUCTIONS  FOR  USE  OF  “MOMENT” 

The  following  are  the  instructions  for  use  of  the  moment  program. 

(1)  Arrays  of  size  up  to  100x16  may  be  analyzed  with  the  card  input 
version  of  the  program;  up  to  100x100  with  the  tape  input  version. 

(2)  Arrays  for  the  card  input  version  are  punched  onto  cards,  one  row  per 
card,  as  five-place  floating-point  numbers. 

(3)  As  many  separate  arrays  as  desired  may  be  analyzed  in  on*  computer 
run.  Each  may  be  solved  once,  or,  if  desired,  any  number  of  times,  with 
the  starting  ordering  chosen  at  random.  In  the  latter  case,  the  overall 
averaged  solution  is  given  as  well  as  each  individual  solution. 

(4)  The  following  data  cards  are  necessary  for  the  program: 

(a)  A  card  with  8  random  integer  digits  in  columns  1-8.  This  is 
always  the  first  card  of  the  input  data  deck,  and  is  required  to 
initialize  the  random  number  generator. 

t.  The  only  oodUtatioo*  neceuuy  ue  to  the  tota  of  the  Input  del*. 
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FIGURE  E-5 .  Flow  Chart  for  Moment  Ordering  Algorithm 


110 


COL 


l?!??,,l’W-i ’Sovrw^'t ’ST s 3 ?.nrv > 


(b)  A  card  with  “FINIS”  in  columns  1-5  as  the  last  card  of  the  data 
deck. 

(c)  In  between  the  previous  two  cards,  separate  packets  of 
cards,  one  for  each  array.  If  an  array  is  to  be  analyzed  m  the 
random  mode,  the  first  card  of  that  packet  must  contain  the 
number  of  random  tries  (all  integer)  to  be  carrk-I  out  in 
columns  1-4,  and  the  work  “RANDOM”  in  columns  9-14.  If 
this  option  is  not  to  be  exercised,  this  card  is  merely  omitted. 
The  next  card  (therefore  the  first  card  for  the  one-time-only 
option)  contains  the  number  of  rows  (integer)  in  columns  1-4, 
the  number  of  columns  (integer)  in  columns  5-8,  and  the  name 
of  the  array  in  columns  9-80.  This  is  followed  by  the  cards 
containing  the  array  proper. 

(d)  As  many  packets  of  cards  for  individual  arrays  as  desired  may 
follow  each  other.  A  sample  data  deck  for  the  card  input 
version  is  shown  in  Fig.  E-2.  (A  deck  for  the  tape  input  version 
could  be  identical  except  that  it  would  not  have  the  cards  with 
the  arrays  themselves  punched.) 


I 

82340795 

3  5SAMPLE  DATA 

1.0  3.0  0.0  0.54  11.3 

8.4  2.7  6.11  12.04  18.1 

0.0  14.0  111.  6.08  10 

10  RANDOM 

4  2MORE  SAMPLE  DATA 

1.0  2.4 

1.0  3.7 

0.0  2.15 

1.0  6.0 

2  3MORE  SAMPLE  DATA 

4.1  7.9  6.1 

4.2  0.0  7.0 

FINIS 


-♦-RANDOM  NUMBER 
-♦-HEADER  CARD,  ARRAY  I 

|  ARRAY  1 

-4- RAN  DOM  CARD,  ARRAY  2 
-4-HEADER  CARD,  ARRAY  2 

|  ARRAY  2 

-4-HEADER  CARD,  ARRAY  3 
!  ARRAY  3 
-4- LAST  CARD 


FIGURE  E-2.  Sample  Data  Deck  for  Cord  Input  Version  of  Moment 

Ordering  Algorithm 
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C.  “MOMENT*  PROGRAM  OUTPUT 

The  program  output  consists  of  the  following: 

(1)  For  each  array  analyzed,  the  program  prints  out  a  complete  copy  of 
the  array  as  read  in,  numbering  the  rows  and  columns. 

(2)  If  only  a  single  analysis  is  to  be  done,  the  program  first  prints  out  the 
initial  value  of  the  correlation  coefficient,  R,  and  then  after  each  row  or 
column  operation  prints  out  the  entire  array,  numbering  the  rows  and 
columns  appropriately  and  giving  the  new  correlation  coefficient.  When 
it  reaches  a  solution,  it  prints  the  array  and  the  message  “THIS  IS  THE 
SOLUTION.’’ 

(3)  If  the  random-ordering,  multiple-attempt  option  is  being  employed, 
each  time  a  new  solution  is  found  the  program  prints  out  the  order  of 
the  rows  and  the  value  of  the  correlation  coefficient.  When  it  has 
finished  the  appropriate  number  of  attempts,  it  prints  out  a  complete 
copy  of  the  most-commonly-found  solution  and  lists  all  of  the  solutions 
found,  giving  the  correlation  coefficient,  the  number  of  times  found, 
and  the  order  of  the  rows  for  each.  It  then  lists  the  overall  averaged 
solution,  giving  the  ordering  of  the  rows  and  columns,  and  the  average 
position  (with  the  RMS  deviation)  of  each  row  and  column. 


If  the  program  encounters  an  unstable  array  structure  (i.e.,  one  in 
which  no  solution  may  be  found,  but  rather  a  cycle  of  states  occurs 
which  will  repeat  itself  indefinitely),  it  takes  that  as  a  solution,  but  first 
prints  the  message  "THE  FOLLOWING  SOLUTION  IS  NOT 
STABLE.” 
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MOMENT  ORDERING  ALGORfrHM  PROGRAM 

PROQRA*  hqm£nT 

trace  arrays 

STEPM£n  OfUTsCN.EXT. 358, SCO 
COMMON  AllOO.lj) 

DIMENSION  NASNAY(SOOO) 

NSOUS«<3»Hl 

uSOLSv2j,iOO) 

M(lO0),Na*),MT(ini)),AT(l6)iMTO(100i,HTO(i6) 
nHjiooj.oKOjioojjONtiii.nNOiieiiiEXTfR) 
ROlD(l0f>)>NOtD<10Q] 

RSAVCdOO) 

AND  N  TO  RQnS 


COLUMNS, 


DIMENSION 
DIMENSION 
DIMENSION 
DIMENSION 
DIMENSION 
dimension 
I  AND  m  REFER 
CPS«l,p.A 
TLEFThq. 

Tn0hicL0CKTm«TLEFT) 
lFnLEFT.LT.t600,  >  SO  TO  14 
CALL  Slu»P(NaHBAY, 5000,0) 

CONTINUE 
R!*D  %9 \ NQUM 

f  opm*t  « 10 ) 

RANBflANIFJNOL'M) 

CALL  9FA0JN(MC,MSlZE,NSm,IEXT,M,ig,NHC»NUM,NMAX> 

TF*  MC.Fo.l )  CALL  »ANORr(M,MSl26jNHC) 

R"«-2. 

NL2ST«o 

CALL  CC(MSIZE.NSI2E,  «,N,R) 

IF(  MC.  K'C.l  ) PRINT  10e,H 
CONTINUE 

Call  Ntn<  M,N,MSI2E,NSIZE,»,MT0,NT0,DN,DN0.MT.NT,IEXT.HC| 
CALL  DtJHPJ  S) 
iMMC.NE.itPRlNT  lOg»R 

Call  mIoj  m,N,mS1ZC«NSIZC«9,MT9iNT0,Dm<Cm0,mt,NT,IIXT.mC) 
call  dijmpis) 

XFtMC.NE.UPRlNT  108, R 
FO»M*Tt,/«  R«*#rl0.5) 

XF«NLlST,CO.«t  GO  To  9 
DO  11  lil, NLIST 

IPI  ABSlR-RSAVLd)  )  .GE.CPS)  ao  TO  it 
IF< I. NP. NLIST)  PRINT  112 

format («  the  following  solution  is  not  stable*) 
go  To  11" 


MaN,hT,NT,IE*T  > 
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n  Old  <  num ) «i 
15  00  1#  T.J.MSIZE 
is  MsoLSifcuH,n*M( j» 

DO  IS  I.l.fcSTZf 

is  nsoLS(Mii*,n«Nix» 

5  PRINT  s;jHT(n,I«l,HSIlE> 

3  ^ORM*T( ,1 x,3sxs/ix,35X3/ix,  JO  in 

print  4 ; r 

<  Format <4h  P»  »rlO,5) 

21  ir (NMC.ttE .NHAXJ  GO  YO  a 
GO  TO  90 

t  KOLD«X»iNOLD<X»*l 
GO  TO  9 1 

8  CALL  $UM( XfXT.MJ*, "SIZE, *Sl£E. OHO, n*0,MSaLS»NSOLS,«Ol.D*iilOLO, KMC 
GO  TO  1 

EnO 

SU8«0UTx\r  REAOINtMc.HSIZE.KSIZE.IZXT.M.K.igMC.KUM.KKAX) 

COMyON  A(tOO,l*j 

DIMENSION  M<Tft!'>,K«l4)»MTti00j,|kT<l4),itXT<9) 

DaTa<I0aN»4MRANDOK> 

«C"0 

5  BE*C  a,MSXZtr,K5IZC,IEX* 

2  F0aMAT<2l4,.)A8J 

K*»C9° 

KOM»0 

IKIBA^.EQ.IEKTII)  J  GO  TO  1 

DO  4  I.l.MSIZE 

«|AD  5,f AJI.J>.J»1»KSXZE» 

«■ CONTINUE 
5  romATtxArs.?) 

DO  «  I«1.NSIZC 
t  KJU*-! 

do  7  m.Mjxze 

7  *IU«X 

CALL  PPINTACMSIZE.KSIZE.  M7N»M?,NT.I|xn 
PHINT  1 1  *  If  XT 

11  f ORMAT(//iflMl5CLUTXOKS/lX»t40/l 

MjTURN 

1  NHAX.«5IZ| 

MC«1 

GO  TO  J 

end 

SU9R0UT1NE  RaNO«D<M,MSXZE«NMC) 

COMMON  a(10041|) 

DXMEKStOK  100, 

DIMEKIIOK  NANH10O) 

9  DO  ir  T.1.M5IZC 
Bans  tn.RANirt-l) 

10  CONTINUE 

Call  SORTstBANS.M.Msm.i) 

NMC«NMe«l 

fitTUBN 

end 

5UIR0UTXNI  SUM! XEMT,NUM.*t«IlC.NlIl|. 090. ONOAMtOLJ#NfOL5, MOLD. 
lNOLD.NMCl 
COMMON  41100,191 
OXMCNIIQN  mx 1 1 00 ) 


H4 


--37  '  ftvWiUsSjH^S’^.  -  ViO  ;f'.v--'^T-^v 


OIHENIION  **SOLS(2j,iOOj 
0IMENSJ0N  NS0LS<2S,H) 

DIMENSION  *VC(iOO),jj£V(10fl  t  f^OROERli00) 

DIMENSION  BOLDfiOO),NOLD(10*l 

DI^ENSIO*'  TE*T<9),ht<100j4nT<14)#dhOCmsI2E)iONOCNS2ZE) 
DIMENSION  NxdA) 

12  IS0L*1 

If  C  NUN.I.T  .  2 )  GO  TO  26 
DO  22  T.2.HUM 

22  IPCN0LP«I  I  .GT  .NOLDHSOL )  >  IS0l«I 

26  Hi«fl 

.  DO  33  Til, “SIZE 
33  Hx<i>.HSOLStlSCL#I) 

DO  34  iil.NSXZE 
3*  NX C I >«NSOLS< ISOL, I> 

CALL  P0INTa<mSIZ£,NSI2E,mx,NX,mT,NT,I£XTI 
print  *? 

35  -*ORmaT(///i  THIS  IS  The  “OST  COHHON  SOLUTION*} 

PRINT  13,11X1 

13  f ORH*T <1h1 /•  SUMH*«v*/1 X,6A6^/«X,*B*,6X,*TIMES  fOUNU 
1/J 

,  DO  25  III , N5I2E 
25  ir  X DHO  t  T » .|(3.  o  , }  N2iM2*l 
mlOim2*«mjj:?-»i2I/J 
mhI«h5XTE*1-<H$I2e-m2)/3 

HG*MZ*1 


“OH  OROEH*/l 


Nh2«0 

,  DO  H  III ,NS12f 

14  ir<DNO«n  . C j .  0 , :  NH2«nmz*1 
Do  23  IIil.NUH 

DC  24  I.i.MSTZE 
NX<I»»MSOL9(TX,IJ 
24 

call  iobtichx.nt^size,*!) 

VQtciO 

DO  2T  IiM*,mlO 
HVAL*HSOLS( isol  » NTC I ) ) 

IrCHVAl.LI.HlOj  H0K.M0«*1 

2?  Ir<H«'iCf0l.MMXi  hoh*hoh«i 
DO  2i  I|HHI.MJI2C 

HVALMSOLICISQL.HTC  in 
ircnvAi .li.hloi  hcxihox.i 

3e  ircHVAL.st.HHXj  hok*h<3<„:. 

I»CHOK.OT.O)  GO  TO  Jl 
_  DO  2*  IiJ.msixc 

29  tr«HSOLSCXI.I».QT.Ht>  HiOLS<II,2>«HSIl|*H|*i-HSOLS(II,I) 

,4  DO  3t  Til.NSHE 

3?  It  C  NSOLSf SCI ,T  I  .P. T ,  N«Z )  NSOL»< X X *1 > •NSXZE*kNI*l.NSOL$C IX . I J 
DO  JO  TiJ.HSm 
hX<I»«hR0L9CIX,X» 

10  HTCIlHt 

CALL  IORTtCHX,HT,HSX2E,*tt 

11  PRINT  »?,ROLOC!IJ.NOL5CII)#CHTCJ»,J«l,HS:fE) 

12  roRHtTclX,rl0.5,X«,5X,S2Xl/tax,l4x3/l«x,34i3) 

23  CONTINUE 

print  jj 

15  EORNATCA/r*  AVERAGED  SOLUTION*///*  Squ  AVERAGE  DEVIATION* 

1/1 

AJINHC 

DO  16  !*l,x$X2E 


^SattKsaBKaeiBw  _ _ 


AVE ( I )»0  , 

DO  1 6  J.1.NUH 
A**MS0LS< J.I> 

AV*N0L0(Ji 

AVE<I|«aVE»H*AX*AV/AZ 
*6  N0RdE»(D«I 
DO  1?  I.I.mSTIE 
DCV(U«ft. 

do  la  j«i,num 
AX»MSOL'S«J»I> 

Ay»N0L0  -t  J » 

la  0EV<I)«1)EV<I)*AY*<  AX«AVEf II  >**2/AZ 
17  0EV(I>B!t5S?(nEV(in 

CALL  S0»Yt (AVE, HORDE". **SI«i-i> 

DO  lal.KSTZE 

"BINT  30,N«RDER<U»AV£(I>,oev(NO»DCR?l>  > 

20  PoRWAT(ia,?FtO ,2» 

*9  continue 

PRINT  M 

38  f0RmaT( ///•  AVERASED  SOLUTION*///*  COL.  AVERAGE  DEVIATION* 
1/1 

DO  34  I»1,NSIZE 
A  vE  ( I )  »0  . 

DO  3«  Jal.NUM 
AX*NSOlS(J.I) 

AV*NOLfl( Jl 

AVE<I).AVE(I)*AX*aV/AZ 
36  NORD£R( I ) «t 
DO  41  i;i,NSIZE 
DEV(I).0, 

DO  39  Jal.NUM 
AX»NSOLS(J.II 
A V*NOLP ( J I 

39  Oevm«DfV(I>*AV*UX-AV£CU>**2/AZ 
41  D(Vm»90RT»0tvmi 

call  sortk ave.norde».nsiic»-1) 

DO  40  T.l.NSIZC 

PRINT  ?0,NORQER(I)*AVC(I).D|V(NORDCR<n  i 

CONTINUE 

PRINT  Ji 

21  EOPHAT(iMU 
END 

SUBROUTINE  PRXNTA(NSI2E.N»m.  **. 4. NT.ItT. X£XT > 

COMMON  a ( 1 0  0 , 1 4  | 

DIMENSION  MXdOO) 

DIMENSION  »*<»*SX2C>***taiSl*Ei«»AT<A»SIZC*«Nt«N*IZC» 

DIM£NS’^#^  it* T(9l 
DO  3  1*1. NSXZI 
MXIXIUMdl 

3  "Tlilil 

CALL  SORT||MX.MT.NSlZC.*ii 

do  s  iii.wsm 
MXII)*N|II 
S 

CALL  lORTf(HX.NT|NSIIC»*H 
PRINT  18. If  XT. NT 

10  roRHAT(lMi.9A«///lN  .6X.1RIT/) 

DO  10t  Iti.MSlfE 

PRINT  teO.MTtt  I.  lA<MT(I).NTTjn;j«t»NIXlll 


too  CONTINUE 

200  FO»MAT<iH0,I*,l4r7.3) 

END 

SUBROUTINE  HID t  H,N,MSIZ£,NSIZI#R,mTO,NTO,DM,DNO,NT,nT,HXT,MC) 
COMMON  A«100,i»J 
DIMENSION  mx«100) 

DIMENSION  M(MSIZE),N(NSlZE»tMT(KSIZE)«NT«NSIZEI 

DIMENSION  mTO(msIZE>*NTOiNSIZE>TDM(mS!ZE>iDMO<mSIZC) 

DIMENSION  ITABLE<H>,ITir<100),NfSAVE«l°0|  ,IEXT<*> 

ERS«1  ,f«5 
DO  3  X«1 #  MSIZE 
MX ( I ) «M( I ) 

3  mT(I)«i 

CALL  S08T2(MX,MT,MSIZE»-U 
DO  5  1*1 , NSIZE 
HX<I  J.Nm 
5  NT(I)«I 

CALL  SORTED, NT, NSIZE, -1> 

DOllO  f.l, MSIZE 
110  MTOtlUMTllJ 
D0115  J.l. NSIZE 
US  NTO(J»«nT(U> 

DO  20  Iil, MSIZE 
Dl«0 , 

*1*0. 

DO  15  J*1 » NSIZE 
Ol«Dl*A(MTm»NT< J>  } 

XN*J 

Xl«Xl*XN*A«MT(U,NT(jJ  I 

is  CONTINUE 
0«<I>*V1/Dl 

20  CONTINUE 

,  DO  12  1.1, MSIZE 
12  D“0«MUm«DM<U 

CALL  OBOERlDM. MSIZE, MT.H.mT#} 

9iE5T«»2, 0 
««*MSI7E-1 
DO  21  I «1  ,KK 

HIE  ( 1 1  «o 
J*X*1 

xr(DMO(NT(xn.co«oNo(**T(umiTii(Z)«i 
xr(DMO(MTm».eoto.)  xtieii>»o 

21  CONTINUE 
ITIC  *  msxzc  »«0 
NTIED*i 

1*0 

22  1*1*1 

XHXTXE|X>.Nt.0>00  TO  83 
Xf(NTXPB.N|,i>00  TO  2« 

XVtX.ftC.NNtao  TO  2S 
00  TO  22 

23  ntieo*ntiio*i 

00  TO  29 

24  ZTASCCf 1 

XMNTXfO.IT.*>  MINT  94 
Xf ( NTXfo  .  ST  ,  * )  NTIED** 

14  roRMAT(.MQiC  TMAN  «  RONS  01  COLUMNS  txid  WHILE  doing  ALOOSITHW,  f 
lIRST  «  IN  TXf  WfRf  RCO*OE«fO»  OTHE»*  LC^t  XN  OL0  0»0E**> 

DO  2*  II»i«HSXZC 
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I 

I 


2?  MfOdn.MYcin 
B8£STi»2.0 
nfact»i 

DO  28  TI»1.NTXED 
2j  NrACT*MF*cr*ix 
DO  29  H«1,NFACT 
CALL  PCRMUTE(NTIED*ITAflLE> 

DO  5°  JJ*1.NTIED 

50  MTO<I»MTI€D*JJ (■HT(I*NTXED*JTABLE«JJI  I 
DO  31  JJ«{ jMSXZE 
MX< JJ)«MTO<JJ J 

31 

CALL  SOflT2<HX,M,MSIXE,-l) 

CALI  CC  t  MSI2E  *  NSIZE .  m.n.N) 
iPtR.LT  .<R&EST*£PS>)GO  to  2f 
B8EST«0 

DO  32  JJ»1.MSIZE 

32  mt$aV|(JJI«mTONJ) 

29  CONTINUE 

DO  33  JJ»1*NSIZE 
MT( JJ)*mTSAVE( JJ> 

MX<  JJ)«mT»*VC< JJ ) 

33  MNJ)«JJ 

CALL  S0RT3(MX«m*mSI2E»-1) 

NTIED*1 

ird.LT.xioao  to  22 

25  R«BBEST 

IFt  B .EO  .*2  •  0 (CALL  CC»NSI2E«1$I*t*  M«N#R> 

IF(MC.F0.1>BCTIIRN 

PRINT  10.IEXT.NT 

10  roRMxT(iMl,9A8>/^lH  ,*X,1*IV» 

DO  100  i»i,H5l2E  .  .. 

print  joo.mtc: j, (a(mt«i) .ntu> >?d«i#*ix2i»#DM0{NT<x)> 
100  CONTINUE 

200  F0BM*T(iM0,l2#l«F7.3iF3.2> 
return 
End 


subroutine  nio<  m.n#nixzi«oisizi»b»mTo»^t0.dn.dno,**t,nt,xcxt,mc) 
common  A< 100,1*) 

DIMENSION  MP| 1*0 ) 

DIMENSION  M(MSII|)«N(N3lIEI|MTtli8IZl)»NTtNSXlEI 

DIMENSION  MTO(MSII£).NTO(N«|2E)IDN«NSllE),DNO«NSXIE) 

DIMENSION  IT able <  »XTIE<1*®> .NTS  A VE ( lOO ) ,IEXT < • ) 

EPStl.f.5 
DO  3  lal.MSIZC 
MXII)bm(I) 

3  MTdltl 

CALL  l08Tt«MX.MT,NSllE.*l) 

DO  5  I ■  1  • NS12C 
MXlI)aN(I) 


5  NTII)»I 

CALL  16RT2(MX.NT|NSI2E««1I 
DOllO  I.1.MSIZE 
110  MT0(I»«NTII) 

colls  J.l.NSIIl 
Us  NTO  ( J  )  «NT  (  J  ) 

DO  20  J.i.NSIH 

Cl"®. 

*i«o, 


I 

I 

l 

I 

I 

I 

I 
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DO  15  I.1.MSI2E 
D*«DU*|MT<I)  »nT»  J) ) 

XN»I 

xi«xi*xn**(mt<i)«nt(  j>) 

15  CONTINUE 
ON< J)«*1/D1 

20  CONTINUE 

DO  12  Ial 4 NSIZE 
12  OnO(NT(i>)»on<I) 

CALI  OBDrB«DN»NSIZE,NT,N,NT9) 

«8EST«-2.0 

XA»NSlZE*i 

DO  21  !.1,«K 

ITIE ( I >  »0 

J«l*l 

IF(DNO(NT(I> j.EQ,ONO«NT( J)» ) 1TIEC  X ) «1 
If(DNO(NT<I)>.EO.O.)  ITIE ( I > »0 

21  continue 

ITIE (NSIZE >»0 

NTIED»1 

laO 


22  I«In 

ir<ITlE<IJ .NE.OJSO  TO  33 
ir«NTlE0.NC.i)80  TO  24 
irCl.GE.KKiaO  TO  25 
GO  To  32 

23  NTIEO»n jlECn 
GO  TO  ?2 

24  ITABLE<1>«-1 
lrfNTlED.8T.fl)  PRINT  34 

ir(NTieo.GT.fl)  ntied««  . 

34  EOBM*T(«no»E  Tn*N  9  ROWS  01  COLUMNS  TIED  WHILE  DOING  ALGORITHM, 
1IRST  9  IN  TIE  «E»E  REOROEBCOi  otwebs  leet  IN  old  0»DER*) 

DO  27  II»1 1 NSI 7£ 

27  NTO(II).NTIII) 

R8ESTa«2, 0 
NfACTii 

DO  28  II-l.NTZED 
2e  NrACT«NFAOT*II 
DO  2fl  II»l,Nr*CT 
CALL  PERMUTEfNTXED.ITAILC) 

DO  JO  Jjal.NTltD 

10  N70(I*HTIfD*JJ)«NT(I*NTIEO*lTABLI< JJ) ) 

DO  31  JJ*1, NSIZE 
MX( JJ ) »NTO ( JJ ) 

31  N  |  J  J  )  ■  J  J 

call  sort2(mx»n«nsize(>i) 

Call  CC(M«IIE*NSXZE.  N.N.ll 
Ir«B.LT.«BSESl*EE*>>GO  TO  2f 
BjESTaB 

DO  32  JJ>1>NSX7E 
12  NTS*V|<JJ)*NTOIJw') 

29  CONTINUE 

DO  33  JjaitNSXZE 
NT(  JJI.NTUVf  ( JJ) 

HX( Jd)»NTlAVE< JJ) 

33  N( JJ)«JJ 

CALL  ipRT*IHX»N,N|IZE.»l> 

NTXEOtl 

ir< i.lt.wkigo  to  22 


JW  <  '  •  4V-? 1  " 


jr 
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25  R«RBEST 

IP:s.£B.-2.0)C*LI  CC<NSIZS,*SIZLi  “*«.«> 

IHMc.  CO.  IIRCTURN 
4.  PRINT  10.ICXT.M  A 
10  F0BM*T<lHl,^*8///iH  ,2x,ls!?/) 

DO  1D«  T»i,MSl2E 

PRINT  20('<MT(IJ<(*(RT(I)4NTfJ)17J«i»R5l2EI 
180  continue 

200  fOB*AT<iMo,I2#l«P?.J> 

PRINT  200tXI«(ONO(NT(<jH  ,J»1»NSIZE  > 

return 

END 

subroutine  order <d*msize»mt/m,mtoi 
COMMON  AilOO.l*) 

DIMENSION  MX(100) 

DIMENSION  0<MStZE)  »hUmsIZE>#M(HSIZE),mTO«  MSI2E) 

Caul  sorT2<d.mto1msize,-1i 

DO  60  I.I.MSIZE 
MT(I)«RTO(I) 

60  M(l).i 

CAUL  S0RT2JmT0,M,MSIZE»-U 
End 

subroutine  cce mSize.nsizEi  m,njR» 

COMMON  a<1°°.16> 

dimension  M(MSIZE*<N(NSIZEI 

dimension  XLIST(1O0),vlIST(*00j 
SX«0  , 

sv»0. 

ssx-o. 

SST-0, 

SPXV«0. 

SR"  0  i 

Xum«hSIZE 

xunansiie 

XUM«1,/XLM 

xuN«i,mN 
DO  1  J»1 # NSIZE 
X»N( J) 

1  VuIST( J)»X*XUN 
DO  2  Xal.MSIZC 
XlM(I) 

2  *UI8 T 1 1  j«>!*XLM 
DO  100  j.i.MSIZC 

z»xuisr<n 

DO  100  j*I iNSIZE 
ZIMYLIST(J) 

TR"P*ZZ 
SX"SX*XP 
Sy»SYayp 
S|X*SIX«IMI 
SSV»SSv*YP*ZZ 
SRXVi3PXT*XP*ZZ 
100  SMSP*P 

SIDXa«SX»<SX»**2/|P 

S8DT«lSy*<SV»«*2PIP 


SfDXY«9fXV»SX«SY/5f 
H»»SQ«TfSSBX*SSO'M 
If  enR.ro. o.  >  r«o.« 

If (RR.K'E.O.)  H«8fDXY/RB 
EnDSENOt 
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THE  MOMENT  COMPRESSION  COMPUTER  PROGRAM 
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THE  MOMENT  COMPRESSION  COMPUTER  PROGRAM 


A.  INTRODUCTION 

This  Appendix  describes  the  computer  program  MOMCOMP  which  implements  the 
gradient  algorithm  in  Chapter  Ill  of  Part  II.  The  main  program  logic  is  shown  in  the  flow  chart 
in  Fig.  F-l. 

B.  DATA  RESTRICTIONS 

The  program  can  handle  an  indefinite  number  of  non-negative  matrices,  each  up  to 
75x75  in  size.  The  data  packages  for  successive  matrices  are  stacked  one  after  another. 

C.  FORMAT  FOR  INPUT  DATA 

The  listing  of  MOMCOMP  includes,  for  illustrative  purposes,  the  data  package  required 
as  input  for  the  second  example  in  Chapter  III  of  Part  II.  This  is  a  10x10  matrix  and  three 
starting  points  are  requested  for  the  row  or  column  optimizations. 

The  data  input  package  for  each  matrix  consists  of  two  types  of  cards.  The  first  card 
type  contains  the  three  variables  N,  M,  IRAN  punched  according  to  format  (315). 

Here 

N  =  number  of  matrix  rows 

M  =  number  of  matrix  columns 

IRAN  ®  number  of  randomized  starts  for  the  optimization  hi  each 
direction. (Row  and  column  optimizations  therefore  consume  2  ■ 
IRAN  starts.) 

The  second  type  of  data  card  is  used  to  read  in  a  single  row  of  the  input  matrix.  Successive 
rows  are  punched  on  distinct  cards,  with  each  card  employing  format  (601 1 X 
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If  the  matrix  contains  up  to  60  columns  exactly,  N+l  data  cards  will  be  required, 
namely  one  of  type  1  and  N  of  type  2. 

D.  SUBROUTINE  DESCRIPTION 

The  portions  of  MOMCOMP  include: 

( 1 )  Main  Program,  which  controls  execution. 

(2)  DATA  subroutine,  which  reads  input  data. 

(3)  PERMGEN  subroutine,  which  generates  randomized  permutations  for 
use  as  starting  points. 

(4)  ZMIN  subroutine,  which  minimizes  the  sum  of  either  the  row  or 
column  moments,  using  the  iterative  gradient  algorithm  described  in 
Chapter  111  of  Part  II. 

(5)  LAP  subroutine,  used  by  ZMIN,  which  solves  linear  assignment 
problems. 

(6)  Z  function,  which  computes  the  objective  function,  taken  as 

for  optimization  of  column  order 


for  optimization  of  row  order 
Z  is  the  root-mean-square  moment  arm. 

(7)  Subroutines  INITCOL  and  INITROW,  which  prepare  the  dat)  needed 
by  ZMIN  for  minimizing  Z.  The  required  data  are  the  W‘s  and  H's  in  the 
expression  (see  Eq.  8-9)  in  Chapter  III  of  Part  II). 


where  L  =  N(M)  for  optimization  of  the  column  (row)  order,  and  where 
jr  i ,  a^, . . .  ,ttl  denotes  a  permutation  of  the  L  columns  (rows). 

E.  PROGRAM  FLOW  CHART 


Figure  F-l  is  the  computer  program  flow  chart  for  the  Moment  Compression 
Algorithm. 


FIGURE  F-l .  Flow  Chart  for  Moment  Compression  Algorithm 


F.  ERROR  MESSAGES 

(1)  The  only  error  message  from  ZM1N  is  that  convergence  has  not  oc¬ 
curred  after  100  iterations  (each  iteration  is  one  gradient  step  and 
involves  one  linear  assignment  problem). 
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(2)  Five  possible  ERROR  messages  may  occur  in  LAP: 

Type  1 :  Step  1  has  been  unsuccessful  at  covering  all  zeros 
Type  2:  Same 

Type  3:  Matrix  element  for  the  linear  assignment  problem  is  negative 
(should  never  occur) 

Type  4:  Step  2  fails  to  find  a  primed  zero  in  the  indicated  row  (should 
never  occur) 

Type  5:  Same  as  Type  3  error  (should  never  occur). 

G.  TIMING 

The  following  times  should  be  multiplied  by  twice  IRAN. 

(1)  10x10  problem  consumed  roughly  3  seconds  for  optimization  of  row 
.  (or  column)  order  for  each  starting  point. 

(2)  16x16  problem  consumed  roughly  10  seconds  for  optimization  of  row 
(or  column)  order  for  each  starting  point. 

H.  COMPUTER  PRINTOUT 

The  computer  printout  from  MOMCOMP  consists  of  the  following: 

( 1 )  The  input  data 

(2)  The  number  of  the  starting  point  (ranging  from  1  through  IRAN)  and 
whether  the  row  or  column  order  is  being  optimized 

(3)  For  each  starting  point,  the  sequence  of  permutations  and  Zs  generated 
by  successive  iterations  of  the  subroutine  ZMIN. 

There  is  no  attempt  to  choose  among  the  several  solutions  obtained  by  varying  the  starting 
points  and  no  printout  of  “the”  final  matrix  since  this  is  generally  non-unique.  The  us«r  must 
extract,  from  the  printout,  the  row  and  column  permutations  which  minimize  iheir  respective 
Zs. 

The  computer  output  must  be  interpreted  as  follows.  U  the  row  permutation  is  printed 

as  ir(l) . ir(N),  then  the  optimal  rearrangement  has,  as  its  ir(i)^1  row,  the  row  of  the 

input  matrix.  The  user  is  reminded  that  the  permutations  (?r(l) . sr(N> j  and 

(N+l-*(  1 ), . . .  ,N+  l-t(N)  j  are  equivalent,  each  being  the  other  in  reversed  order. 
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‘  r**iTI*&v,FZw8St!P*  v?- 


*/£-;■  rriVY'  F«;TrvfrF'7ff»?1. ’r^gfr^T^'  yp  < 


MOMENT  COMPRESSION  ALGORITHM  PROGRAM 


io  10  3 

PROGRAM  MQNCOMP 
C  P.SCHWFITZPR  X?i4 

common  w(  100 » .*<75*751  ,D<  75/75!  »C<  75.751 .1«AR.Ni)' 
DIMENSION  7PfH«<lOO) 

C  INITIALIZE  RANDOM  NUHBf  B  GENERATOR 

x»  RAN" fj(J, 1416) 

C  FINDS  IRAN  LOCAL  min  fob  HpW  and  column  PERMUTATION*. 

1  CONTIM'E 

CALL  TTMFM9.17MREGIN  NE"  MATRIX.) 

CALL  Data 
C  CE*DS  P AtA 

Call  TTMF( 39.11MUATA  IS  IN.  ) 

call  initcol 

Call  ttmfmr.ITMINITCOi.  nONE.) 

C  C0«PuTrs  w  and  m  for  column  pebmuTaTIOns, 

Do  2n  T«i .Iran 

1»m 


50 


30 


CALL  PAGfS*<P 
PRINT  50.1 

CaLL  Pfrm5EN( IPERm.l ) 

call  timfjs^.iihpermgen  dons.) 

GENERATES  PERmiTaTION  PF  1-m.STORES  In  IP£Mm. 

Call  TImfi !9,41m;mIn  STARTS  optimization  of  COLUMN  URUE".) 
call  Z“i» (l.ipcrm) 

2»IN  FINPS  optimal  COLUMN  PeRMUTATTON,  STaBTInG  ‘hOm  IPE*m. 
Call  TI«f  t  Sq.SrmZmIn  E'DS  optimization  of  pOLijmN  ouue  r  .  > 
CONTINUE 
CALL  INJTRO* 

call  timf(49.ishinitro»  done.) 

COMPUTfs  w  and  M  for  bow  PERMUTATION 
do  3n  T.1.TBAN 
L«N 

CALL  PAQESKP 
PRINT  50.1 

FORMAT(.PE»GMEN  CALLED  FOR  TIME*# 1 5 1 

call  PrRM8EN( i»erm.l i 

CALL  TTmE( S9,ilHPERMaEN  DONE.) 

CALL  TiwF(S9,3(|H?MlN  STARTS  OPTIMIZATION  OF  »0«  OwDtR.) 

call  iuiN(LitPfRMi 

CALL  TimF(S9,J4WZMIN  ENDS  OPTIMIZATION  OF  bow  ORDER.) 
CONTINUE 


isFjs£^ 


i  be^. 
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tat&hA 


\  1 


EJNDS  OPTI“AL  ®Ok  PERMUTATION.  STARTING  PNOM  IPENh, 

go  to  •■ 

END 


Subroutine  data 

Common  «<100>»m<75#75>.D(75»75»iC{75*75>»JRAN,N»I' 

Call  PASPSKP 
W£*U  loi^.M.lRiN 
EqRmaT(3I5) 

PRINT  15#N,m,IRAN 

FURMATf.NER  CASE.  N.H.JRaN  •  *<4l5) 

DO  1«  Ta1  i  N 

READ  2n;  (0(I.J)iJ»l«H» 

P0BHAT(7Pl°.S! 

OJMENSTO*1  10(100! 

DU  25  I*i.N 

READ  21 ; (I0( J) , J»1#M) 

E0»MAT(6n!l) 

no  23  J«l,“ 
uii.ji*  to ( j i 
CONTINUE 
PRINT  «0»<  « 

f  ORmaTj  y*0»IfiINAL  Matrix#^/*!  J  OtlfUl  */|  215**  1  o  ,  3 » > 
PRINT  40 

FijRmaT(.  DATA  all  READ  in  «) 

return 

END 


Subroutine  initcol 

COMMON  *< too » »M 7J. 7S>.0< 7Sf 7S> *C<75.75> ♦!»**.*»* 
DIMENSION  S(i°0> 

DO  5  I»1 » N 
S|l)»0. 

DO  3  U»1,m 

SIX)  a?iT»*  D(I«J» 

CONTINUE 
DO  15  J.I.m 
«*•  J  >  *0  . 

00  10  lal.N 

w|J)a  i«  I  J  I  «D  <  I » U  •  /S  ( 1 1 

CONTINUE 
00  30  Jal i M 
DO  1«  * ■ I  *  J 
*•0. 

DO  17  lai." 

M| J»N|ai 
CONTINUE 
00  40  J»1»N 

DO  3S  «aJ.M 
•4 1 J  *  X  |  a 

continue 
Call  PAGEUP 

PRINT  5o,|IJ.X< 

ro*NAT</*fNXTCOl  DONE.  J  *<J» 
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i 

I 


I 

I 

I 


A* 


I 

I 

I 

1 

1 

I 

I 

I 

I 

I 

I 

I 


•§ 
.  $ 

:-$fe 


«5r-U 


3EW^^W5S«n5«Sfl^^KS55W^ >r- 
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C 

C60 


10 

15 


17 

19 

30 


35 

40 

c 

c 

C50 

C 

C60 


10 


PWINT  «  0 ,  (  f  (J,K,M(  J,K»  i  ,K*1  #J*i 
fOHM*T(/«J  K  H  <  J  <  K  )  */ I  215, F 16 . 4 )  I 

«£TUBN 
END 


SUBWOUTJNE  INITROl. 

C0M*0N  u(100J  »H( 75*75) *D(75/75),C(  75*75). X»AN,N#M 
dimension  stioo, 

Do  5  I •  1 ,  H 

S(I)«0. 

DO  3  JbI.m 

S(I)  ■$< t> ♦  n t o* i) 

CONTINUE 

DO  15  J«1,N 

DO  10  T.I.m 

"(J)«  W(J)*D(J*I>/S(I) 

CONTINUE 
00  30  Jal *  N 
HO  lg  4,1, j 
X»°. 

DO  17  I.i.N 

x»x*o(j,n*o(K#u/S(ii»«2 

Hf J,K>«x 
CQNT  IM'E 
DO  40  .)■!, N 

DO  35  *«J,N 

MIJ.K),  WJK.J) 

CONTINUE 

CiLL  PiGFSKP 

pNlNT  50,lfJ»MJ)>,J*l,N> 

F0RM4T(y*INITHnn  DONE.  J  W(J1  */( l5*El|.8> ) 
PRINT  40,(1 <J,N,M( J,KH ,K»1 /N| , J»1,NI 
PORM*T  t  /•  j  4  MIJ.KI  •  / ( 22  5*  F It , 4 )  I 

return 

end 


subroutine  pehm&emipebn*n} 

DI*EN$roN  IPER«(lCO,,xv(loe, 


oo  i®  T,1,N 


*MX)B!>*NDN(0) 

C*u  SOBTJ()M|,IPEBN,N,l) 

return 

end 


function  i<*,i»e«m> 

COPMON  uli00|,N(751T5>,DI75(75>#C(75i7*j)iI»4N,NKO4#HC0L 

DIPENSTON  tPf»P|10®> 

i.C. 

DO  10  jat.M 

2«2*u(u,«ri.o*Tr*Ti>E»M(jn**X 

00  %  «»1*M 
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-  ^  -i- 


5  j«FLOiTF<IPF»»«(  Ji*IP£R«<Kn 

1°  COUTINi'E 

Z«  SQRTF(Z/FlOaTF(H) j 

meTurn 

end 

c 

c 

c 

SUSWOUTl^g  ZMIN(H,IP£Rm) 

common  W(  1  0  0  5  ,w<75,75) ,D(75/7S) /C<  75»7S>.I»AN,nHQ«,mC0L 
DI*ENSTC^IPEPH(100)  ,TcnRPER«100)#KiF»PEB(ll'n)  ,v(lr,>) 

ZNEk»  Z(m.TPc»h) 

PRINT  10,  ( (I,IPE«M(ln  a«i,M) 

10  F0BmaT</*  Z«IN  INITIALIZED  *S  FOLLOWS,  HEKMUTATI0N*/<2lb) I 
PRINT  l5.ZN£w 

Is  FOP*'*!  (*C'flJ£^TIVE  F(J  Na  *,£20.6) 

DO  20  Tt1 ,k 

<0  NERpER(I)«  IP£Rm(I) 

DO  100  LLL»1.100 

C  CALL  TT«P  (1<j,2ahZMIn  8ESINS  ONt  ITERATION,) 

LL«LLL 

C  DO  UP  TO  100  LOOPS  Of  ITERATIONS 

„  DO  21  Tii,M 

21  ICUBPEB(T)*  N£ bPE P ( I ) 

ZOLD*  7NFW 

C  GENERATE  C  MATRIX 

DO  25  j«i,m 
Y|J)«  0, 

DO  25  * *1 1 m 

2i  V(j)»  v<jj.2,*w(J,K)*FLOaTP<1CURPEB(K)) 

25  CONTINUE 

DO  30  J.I.m 
DO  2g  x.t  1  h 

28  C(J.Kj*  *(J)*FL0ATF(K).*2*  U  J)»FL0aTF(K) 

30  CONTINUE 

IL«  LAP(  m.NEwFEB.OPJ) 

IP  (IL.FO  .1)  betupn 

ZnE*>»  7  ( ** »  NeuPEO  ) 

PRINT  aO.lI.ZOLD.ZNC* 

40  F 0Rma T ( #  ZMJN  ITERATION  NgM|ER  *.I3,*,OLD  AND  N£h  OBJECTIVE  FuNcU 

IONS*  *,2e?0 , j ) 

POINT  45.MI,NEwPC»{In,Iil7M) 

*5  FORhaT(«N[w  P|RMUTaTI0N*/(2I5> ) 

Ir  (  Zne*..  70LD  >  100,40,150 
40  uo  70  T,lt* 

IF  t  <  I )  .NC.ICURPCBdl  )  GO  TO  1«0 

70  CONTINUE 

GO  TO  )50 
100  CONTINUE 

150  PRINT  160, LL 

1*0  rO»*ATu ITERATIONS  £*0  at  L90P*JI5> 

«|TUBN 

End 

c 

c 

c 

Twnction  giP)  NiIPEPmot.osj) 

COMMON  w( 100 ),M« 75. 751,0(75175 )iC( 75.75 )  »IRAN,Ntl0B,MC0L 
C  LINEAR  ASSIGNMENT  PROetEM  SOLVER  MITn  n  X  n  c.  c  IS  destroyed. 

C  RETURN  yITM  LA**o  IP  SUCCEMPULF  1  IP  UNSUCCESSFUL, 


c 

c 

cl 

c 


5 

10 


is 


20 

C 

C 

C 

C 


25 

26 
27 


30 

35 

C 

c 

c 


OUTPUT  IS  OPTIMAL  PfPMHTATJO*  IN  I«»EMmoT,  minimum  OBJECTIVE  FC*  IN  OBJ 

PENSION  XPeqcOT»i»i-,#m«oel0O»JjZEQoe500  ># 

1N*HK(  500  ) ,  MgOLCOV(  100  ) ,  mHO*COV<  l!'0  )  I  LIST  I  ?n0« 

OIMENSTON  U(10n)jV<iOO) 

call  timfi s9,iohup  begins) 

Su8«0 , 

ISTEPl.o 

lEfiROR.o 

Sup*  CUMULATIVE  amount  SUBTRACTED  FKOm  a  HOm  or  column. 
pPlNT  i  >,J*i.N)>  1*1, N» 

p0»MATr3T5,El«.7,2l5,El6.7.2l5.  E  i  A  .  7,  ?1 5,  fU.  7  ) 

Subtract  small  eieuent  from  each  rom 

dp  10  i.i.n 


««C(I/  ) 

DO  5  J«2 ,N 

x*  AMI* 1  (X.C'I.J) > 

Sub* 


DO  10  J ■ 1 , N 


C<I.J)f>  C<T,J)-X 

subtract  small  EST  ele-ent  from  Eac*  column 
do  20  J«1 ,N 

X*C( 1,J) 

DO  15  I.?,N 

x*  A M I *.  1  <  V,C<  I,  J)  1 

Sub*  si'b*x 
do  20  1*1, N 

C ( I, J  ) ■  C(T,JI-X 

pRINT  1  i  ( ((I,J.C«X»J> )  ,J*M>  »I*l,N) 

STORE  all  7EROS  In  a  OnE-OI“ENSIONaL  aR*A».  NZEKO  ZtwoS. 

NmTm  ZfnO  TS  AT  IZEPO(M) , jieRO(K) . 

M*RK<*>*  1  IF  THIS  ZERO  IS  STARRED,  »*1  If  IT  IS  PRIMED.  *'J  OThER»ISE 

NZERO*n 

DO  30  1*1, N 

DO  30  j » 1 , n 

IE ( C ( I , j ) )  25,27,30 


pPlNT  ?6,I, j 

E0RMAT(/.ELEMENT*,2l5,*RfS£T  TO  ZERO* ) 

00  TO  *0" 

N/ERO*  nZE»0*1 
IZEROI* ZERO)*  I 
UZERC(*  zrpo j*j 
maRK(NZfRO)*0 
continue 

EORMATl  «IEfiO  PRXNT0UT*/«?»I5)> 

InITIaIIZE  COVERS 

mCOLCOviI)*i  IF  COLUMN  I  is  COV£HEP»u  IE  UNCOVERED.  SIM.  POfl 

initialize  stars 
njtar*a 


c  njtar*  number  or  starred  entries 

C  DO  Initial  STARRINQ  BY  SUCCESSIVELY  CHOOSING  ZEROS  mmicm  h4VE  TmE  LEAST 

C  number  OE  other  ZEROS  in  The  COVERE»EO  PORTION  of  EITHER  ROM  OR  COLUMN. 

C  Then  COVER  That  ROM  OR  COLUMN. 

OO  2000  1*1. N 
mrOmcOv  1 1  )*0 
2000  mcolcOV(T)*0 
2100  ICOVERiO 

M  JNLIN J ■ ?*N 
m:nsum*4.n 

C  AT  END  OE  **SS,  ICOVER.O  h£4NS  ALL  ZEROS  a»E  COVERED.  iCOVfR  POSITIVE 


C  *EaNS  i'Nr0V£REP  ZfROS  fXlSTl  aND  L«£ST  IS  THE  WEST  OP  Them,  HAVING  Th£ 

C  ^  E  WE  S  T  <  mtmoxK'E  )  OTR£R  UNCOVERED  ZEROS  IK  ITS  HOW  0*  COLUMN. 

DO  2200  K«l,K'ZpRO 

C  COMPuTP  NUMflEQ  OP  UNCOVERED  ZEROS  IN  LINE  MlTH  K.tH.EXClUDINa  k  ITSELF. 
Iff mrovcOV»IZE«0»K> > .EC.l  )  30  TO  22 NU 
IP<HCOLCOV<  JZERO(K) )  ,E0.1  )  ao  TO  2?'-'0 
C  K  is  U‘'COVPRFD. 

ICOVE  R«i 
INKROw.o 

INKCOL»0 

C  INKR0W(INKC0L>*N0,  of  UNCOVERED  ZEROS  JNRCk  (  COLLlUN  )  Nl'fH  k 
Oo  21  jn  Hl.NZfRO 

if(izebo<k' .eo.izeRo<li)  ao  to  2120 

ir(JZER0<K>.NE.j7ER0(Ln  SO  To  215'.' 
c  L  IS  IV  Same  COLUMN  as  K  and  L.NE.K,  nQ*  TEST  IF  UNCOVERED . 
ir (MROkcnv(IZERO'L) > .EO.n  ao  TO  2150 
INKC0L«INKC0L*1 
GO  TO  ?t«0 

2120  IF(mC0lcOV( JZ|R0<L) ) .EO.ll  30  TO  2i5u 

irtL.co.M)  ao  to  2150 

c  L  IS  IV  Same  ROM  AS  K,  DISTINCT.  AND  UNCOVERED. 

INKROWrINABOW*1 
215  0  CONTINUE 

iNKLI-'fi  MINOI  INKR0w,iNKC01.> 

InkSum«inkbou*inkcOL 

IF (INKLINI.QT.MINLINE)  GO  to  22 0 (J 

IMINKI.INE.E0. “INLINE. and. I,»*SUK.GE.*INSui'|  GO  TO  2200 
MINLINr.  inklINE 
mINSUm»  INKS'JM 
LBEST.  k 

IF(MinliNE.£O.O.AND.minSUN.LE.I)  GO  TO  2SU0 
2200  CONTINUE 

IFIICOVFR.EQ.O )  GO  TO  2500 
C  UNCOVEOEr  ZERO  AT  LBEST  IS  NON  COVERED. 

2500  MRO*COVfTZERO(LB£ST) )»j 
MC0LC0vj jZFRO(L8EST))b1 
mARk (LPEST).i 

NSTaR«vsTaB»1 

C  PRINT  ?4S0.  IZEROtLBESTl , JZIR0ILBE5T) 

C24S0  FORmaT(*INITULIZaTiON  COVCIS  ZERO  At*.  2j5i 
GO  TO  ?ln0 
2500  CONTINUE 

0  STARRING  0»  ZEROS  DONE .COLUMNS  PROPERLY  COVERED*  UNCOVER  HOnS 
0040  Til .N 
<0  MRONCOVf I >«0 

C  PRINT  S5,MK.IfE«O<K»,JZEH0tK»#NARK(lin*K»l.NIEHOi 

C  PRINT  *0,N5TaR 

Q#0  FORMaT(  •PRELIMINARIES  DONE »  *  »I5.*STaR»,*  ) 

IFINSTar.EO.NJ  GO  TO  400 

c 

c 

c 

C  STEP  1 

100  IST£P1«  IST£P1*1 

C  CALL  TlMfl JS.14MLAP  BEGINS  JT£P1> 

lPIISTERl.aT.SOO)  GO  TO  4 05 

iitepfl.o 

101  XFLAGia 

XITEPFL  •  IST£PFL*1 
IKI3TCM'L.G?.?*NICR0>  go  TO  104 
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C  AT  END  Of  LOOP.iT  180,  IFLAG  aO  IF  NO  UNCOVERED  ZENOS  EXIST,  ;,*D 

c  *i  if  possiblt  uncovered  zc «os  exist,  fxtm*  p*s$es  through  this  loop 

C  APPEAR  T"  PE  NECESSARY  TO  ENSURE  N«  UNCOVEBEP  ?E»OS  EXISU 
DO  1AQ  K*1 , NZ£BO 
KK«K 

!<■  IZEROJK) 

J<«  JZFRC<K> 

Ir(<PROuOOV(IK>*KCOLCOV(jX)».GT.O)  SO  TO  18" 

IFLAG«1 
NARK { K ) b«1 

C  PRINT  125»TK,UX 

Cl25  F0R«aT(  *STE°1  HAS  PRI*-ED  AN  UNCOVtwtL  ZERO  aT*,2I5) 

C  IF  There  no  STaRRED  ZERO  IN  «Ow  IK,  SO  TO  STEP  2.  IF  THERE  IS  A  ZERO 

C  AT  t,  '•OVER  This  ROw  a*D  UnQOVeR  The  cOLUhn  of  l. 

c  repeat  ttll  all  Zeros  are  covered,  then  go  to  step  s. 
no  iso  .NZERO 
IF ( IZE°0  <  L  '  .  NE  .  IK  >  GO  TO  ISO 
IF<L.EC.*>  GO  TO  130 
IF ( “aRk  <L ) . NF  •  T )  GO  TO  130 
GO  TO  15c 
130  CONTINUE 
C  PRINT  1 3s 

0135  FoRM4Tf«?TEP  1  FOUND  NO  STARRED  ZERO  IN  T*IS  RO*,w£NT  TO  STEP  2*1 
GO  TO  ?0A 
150  IpIZERo  t L  ) 

JiJZERO( 1  ) 

C  PRINT  155. I, J 

Cls5  FORhaTj.STFP  1  FOUND  STaRRED  ZERO  in  RO*.  aT  *,  2l5) 

MROwcOv  » ! ) al 

MCOLCOv  y  j  » • 0 
IFLAG«1 

180  continue 

I F ( I F  L  A  S  >  190.190.195 
*9n  CONTINUE 
C190  PRINT  191.ISTEBFL 

Cl9i  FQRMaTi  »STeP  1  DONE  IN*.I5,*  PaSSES.  On  TO  SUP  3,.) 

GO  TO  SO11 
193  CONTINUE 
Ci95  PRIM  19a,ISTEpfl 

ci96  f0rmat(*stepi  "nuone  after*#I5.  «  passes,  resume,*) 

GO  TO  101 

c 

c 

c 

c 

C  STEP  2.  My  KP1«  1*MUNXBES  x 

C  LIST(J)  CONTAINS  The  NUMBER  OF  ZSUBN»1> 

C  AT  210,  STEP?  ASSUMES  7Su8t*Pi>  EXISTS  AND  SEARCHES  FOR  ZSUB<*Pl*l> 

200  Kplal 

C  CALL  TlMFJ39,iftMLAP  BEGINS  STEP?) 

LIST (1 1  a  K* 

210  JTESTa  JZCRO ( LIST ( KPl ) ) 

DO  220  l*1.N?CR0 
IF  I  Mark *L )  .Nc • l>  50  TO  220 
IF(JZEBOIL).CO.JtCST)  GO  TO  2*0 
220  CONTINUE 

C  SEQUENCE  MAS  TERMINATED, 

C  PRINT  j2i,ni,LlSrtI),IZC«OILI3T(X)l,jUH0«LlST(IH»,I»i.KPl> 

C22l  FORMATf  A/al*ERO  SEQUENcr  FOR  STEP  2*/<<l5>| 

C  STaR  PRIMES  IN  SFOUENCE.  UhSTaR  STa*)5  JN  SEQUENCE 
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DO  230  xk.i1.kP1 

IF ( mARk ( L  [ST(NX) )  )  222,230,225 
222  M4PK<tlST'|MK>  )«1 
GO  TO  ?3n 

225  «A*K  <  LIST  J»0 
230  CONTINUE 

C  ERASE  ALL  PRIMES.  UNCOVE 0  ALL  R0“S,  COVE*  EVERT  COLUMN  CONTAINING 
C  A  STARRED  ZERO. 

13 0 24  o  T«1 »  N 
HCOLCOV(T l«0 
240  MRO*COV(IJ»0 
NSTABbP 

DO  250  *«l,NZERO 
IF<maR*<*)  »2a2,250,245 
242  P'aRk  ( k  i 

GO  TO  ?5rt 

245  mCOLCOV(J2ERO(K) I»1 
N$TAfi«  n5TAR«1 
250  continue 

C  PRINT  ?55»NSTAR 

0255  fORMAT<  #5TpP  2  HAS  TERMINATED  kJTn  *,15,  *STaHS*I 
C  PRINT  T5»<(K,I2EH0<K),J?ER0<Kl*MARK!«n»4«i,NZfR0i 

C  PRINT  325,((I#mROkCOV(I),-c9lCQ»<I>>*I»1«N> 

If (NST*b.eO,nJ  GO  TO  4 p 0 
GO  TO  10O 
260  *P1«kP<*1 

LIST(KP1I«L 
ITESTb  IZER0(U 

C  A  PRlMfo  ZERO  IN  ROW  IT£ST  IS  GUARANTEED  To  EXIST, 

DO  27 o  K«t,NZERO 
If  (MARK(K>  ,Nf  ,-U  GO  TO  27° 

If  ( IZEPQ(K) .EQ.ITEST )  GO  TO  280 
270  CONTINUE 

Go  TO  «0i 
280  XP1»  «P1*1 

LIST  t  K*i 1 •  K 
GO  TO  21« 

C 

C 

C 

C  STEP  3 

300  If  (  N'iT AR,  £0  .  N )  GO  TO  400 

c  call  time< s».iahlap  begins  steps* 

C  less  Tm*N  N  INnEPENDENT  STAlRED  ZEROS 

c  CQMPUTr  X,  TME  minimum  UNCOVERED  NUMBER •  IT  IS  STOIUTUT  POSITIVEi 
C  ADD  K  TO  EACH  COVERED  ROW.  SUBTRACT  X  fW0»  EACM  UNCOVERED  COLUMN. 

C  PRINT  R20.NZER0 

C320  FORM aV ( «ST|P  3  BEGINS  hITm*.  15.  *ZEROS.  •) 

0  PRINT  S25,l(IiMROwCOVm.MC<HCO*«X>>«X»i»M 
C325  EORmaTuCOVEB  PRINTOUT*/ < H J* > 

Xa  i,0r*100 
DO  351  J»1.N 

If ( MRQU50  V ( I ) . EQ  «  3)  GO  TO  SSI 
00  350  JbI.N 

If IMCOLCOV(J) .EO.i)  GO  TO  SJ° 

XbAMINI  U.CII.JII 

350  CONTINUE 

351  CONTINUE 

C  PRINT  360.2 

CM  fO«MAT«/.STEP  S  SUBTRACTS*.  C20;8> 
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U 


r 

If 


<? 


370 


375 

3a0 


365 

3<»0 

C 

c 


3900 

39Ci 


C 

C3q0s 

C 


3910 


3V20 

C 

C 


392 

393 


394 

395 
C 

C39« 

C 

c 

c 


IFIX.IF.O,)  PRINT  37° 

FORM*T( //«5TfP  $  ERROR.  UNCOVERED  ELEMENT  NUN-PQSITIVE*  I 

if  (x.le.o.)  an  i o  *02 

DO  3fl0  x»1,N 

IKMOOkcOV(I)  .EQ.0)  go  TO  380 

SuR«  Shb-x 

DO  375  J«1»N 

C(I,J).  C(T> J1*X 

CONTINUE 

DO  3  "  0  J ■ 1 »  N 

IF(MCOLCOV( J) .EO.i)  GO  TO  390 
SuBs  S"8*X 
DO  3B5  I«1*N 

ca.ji  ■  cii»d)*x 
CONTINUE 

DETETE  twE  ZFHPS  WHICH  BECAME  POSITIVE.  These  A«E  PHfClStLV  The 

T.lCE-rovE»EP  ?E«OS 

Kao 

K*K«1 

T.FIK  .at.NZEHO)  GO  TO  3920  , 

lF<<HROueOV<lZeeOi«)UMCouOVUEEftO(*|)).  NE.2)  GO  TO  390: 

PHINT  S9ft5»I2E90tK)*J?FP0|H> 

FO»«aTu  step  X  DELETES  ZERO  at  *.  215) 

DELETE  K-tH  7E°0,  PUTTING  LAST  ZERO  IN  T«IS  SLOT. 

IF(K.EO.NZERO)  GO  TO  Sg  3  0 
m*PK(K).  MaPK(NZEBO) 

h*Hk(N7epot»P 

IZEROtv,.  IZEH0<NZE«0) 

IZE«0<V7i*«0>aO 

J  Z  E  R  O  (  *  1  ■  JZESO ( NZEHO  J 

UZERO<‘'7fR«>*0 

NZEROp^FRO-I 

GO  TO  39^1 

h*Rk ( N7FpO  >  an 

IIERO(N7FPO)aO 
JZERC  ' 7FRO ) a? 

NZEROa  NZfBO.l 

CONTINUE 

ADD  *Ny  new  ZEROS  TO  LIST.  T«ESE  CAN  0NLT  BE  IN  The  UNCOVERED  a*Ea, 
SINCE  all  ZEROS  ON  t.lST  a«E  COVERED. 

DO  395  yai ,  N 

ir(HROWcOV«I» -CO-1 >  GO  TO  3f5 

DO  394  J  a 1 1 N 

TFtHOOUcOV«J>.rO.I)  GO  TO  *9* 

IF ( C ( I » J  )  J  392 1  393(394 
C ( I , J  )an  . 

PRINT  ?6,I.J 
NZERO«N7FRO*l 
HERO  { N7F9O  J  al 
JZERO«NZFRO)a  J 
HaRk (N7ehO>«0 
CONTINUE 
contin'e 

PRINT  194  .NZERO 

FORMAT)  «STEP  3  DONE,  •ill#  *71*08  IN  MATRIX,*) 

PRINT  1  ;  H  II.  J.C<I»J»  .I*1  »N* 

PRINT  15, MK,I7Eh0<5> . JZEB0»K J.XAfittlK}) lKpl .NZfcRO ) 

PRINT  !2«,«  <T,pROnCOV»I)#hC0!.CO9«X>  >.I»iiNl 
GO  TO  TOT 
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DONE  ,  n  STARS  EXIST, 

08J«  Sup 
DO  410  K»t,N 

iPERHOTjtf  j«0 
DO  420  K«l,MZCPO 
ir^APKfK)  .NE.t )  GO  To  420 
IP£Shot  (T?EBO(K  SOfK  » 

CONTINUE 
PRINT  ijn.OHJ 

F0«m4T(.Lap  SUCCESSFUL.  OBJECTIVE  >  *,E20t8//«  QPTIMal  *laCE  r0« 
*1  * ) 

PRINT  44(',t(T.lPP*<"0Tlf)  J,X»1,N) 

FORmaT(2T8) 

P«INT  i  1 1  1 1 ,  J  ,  C  <  I » J  )  ) ,  J«l,  N  )  , I«l, n ) 

UP«° 

RgTURN 


ERROR  “ESSAGES 
IERROR.  TCRRnfl*l 
l£RROR«  TrRROR*l 
IEbROR«  terro«*i 
IERROR*  TERRO»*l 
lERRORt  TCRROR*l 
PRINT  A?n, TERROR 
PORMATf.LAP  rso0R  OF  Type* 
LaP«1 


RETURN 

end 

end 
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DISCUSSION  OF  MEASURES  OF  EFFECTIVENESS 

A.  RELATIONSHIP  BETWEEN  PLEASING  PATTERNS  AND  MEASURES  OF 
EFFECTIVENESS 

The  ultimate  goal  of  any  data-organizing  algorithm  is  the  discovery  of  an  informative 
pattern  of  variable  relationships,  as  evidenced  by  a  pleasing  matrix  appearance. 

Quantitative  measures  of  effectiveness  are  used  as  surrogates  for  pleasing  patterns,  since 
the  latter  concept  is  an  intuitive  one  not  easily  described  in  words.  The  two  MEs  used  in  this 
report,  the  summed  bond  energies  and  the  summed  moments  of  inertia,  were  chosen  with  the 
hope  that  they  would  produce  pleasing  patterns  by  creating  dense  blocks  of  numbers.  No 
doubt  other  MEs  can  be  devised  for  this  purpose;  the  two  proposed  hero  are  useful  because 
they  are  both  (1)  amenable  to  simple  algorithms  for  approximate  optimization  and  (2) 
successful  at  producing  informative  patterns.  Any  other  useful  ME  must  share  these  two 
properties. 

The  algorithms  used  for  optimization  of  these  two  MEs  (the  sequential  selection 
algorithm  for  the  bond  energy  ME,  and  the  gradient  algorithm  for  the  moment  of  inertia  ME) 
are  suboptimal,  that  is,  they  do  not  rigorously  optimize  their  respective  MEs.  Neither 
algorithm  should  be  faulted  for  producing  suboptimal  solutions,  because  the  ultimate  goal  is 
producing  informative  patterns,  not  rigorously  optimizing  the  ME;  the  ME  is  merely  a 
surrogate  for  measuring  the  pleasingness  of  a  pattern.  Indeed,  the  satisfaction  with  the  two 
algorithms  is  based  upon  their  prooucing  data  orderings  which  are  informative. 

It  often  happens  that  several  appealing  data  arrangements  exist,  ah  with  approximately 
the  same  ME  (namely,  near  the  optimum),  and  all  very  similar. 

Consequently,  ties  or  near-ties  among  the  ME  can  only  be  broken  by  a  subjective 
eyeball  judgment  as  to  which  data  arrangement  is  most  pleasing.  Until  the  eyeball  judgment  is 
made,  the  tying  and  near-tying  configurations  must  be  considered  equally  aeeptable.  For 
example,  the  five  solutions  in  Table  2,  or  the  four  solution  matrices  in  Figure  30,  and  the  two 
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solutions  in  Figure  28  must  be  considered  equally  acceptable.  It  is  highly  arbitrary  to  choose 
one  over  the  others  on  the  basis  of  the  numerical  value  of  the  ME. 

In  short,  the  ME  is  useful  only  for  the  first-order  task  of  locating  a  handful  of  good 
arrangements.  The  ME  is  not  useful,  except  in  an  arbitrary  way,  for  the  second-order  task  of 
choosing  among  the  good  (and  nearly  equally  pleasing)  arrangements. 

B.  GENERALIZATION  OF  THE  BOND  ENERGY  ME 


The  bond  energy  ME  can  be  generalized  to  include  bonds  between  matrix  elements 
which  are  not  nearest  neighbors.  For  example,  a  ME  which  weights  the  bonds  according  to  the 
inverse  square  of  the  distance  between  the  matrix  elements  (so-called  gravity  model)  would  be 


ME  =  £ 


£  aii 

ij  rs  ^  y  (*  *  r)2  +  (j  -  s)2 


It  may  be  conjectured  that  such  generalized  MEs,  when  optimized  over  all  row  and 
column  permutations,  are  more  successful  at  producing  tightly  clumped  matrix  elements  than 
the  ME  used  for  the  Bond  Energy  Algorithm,  which  involved  only  nearest  neighbor  bonds. 
There  are  two  objections  to  the  generalized  ME,  however.  One  is  the  significantly  greater 
computational  difficulty  in  optimizing  the  ME  over  ail  row  and  column  permutations.  Once 
the  nearest  neighbor  feature  is  abandoned,  the  sequential  selection  procedure  described  in 
Appendix  D  cannot  be  used.1  Even  more  serious  is  the  fact  that  whet:  diagonal  bonds  are 
included  in  the  ME,  it  is  no  longer  possible  to  optimize  the  ME  in  two  passes,  one  which 
optimizes  the  row  order,  the  other  optimizing  the  column  order.  Instead,  one  would  probably 
have  to  iterate,  as  in  the  moment  ordering  algorithm,  between  row  rearrangements  and  column 
rearrangements. 


Hie  second  objection  to  a  generalized  ME  which  includes  diagonal  bonds  is  that,  for 
sparse  matrices,  optimization  of  the  ME  may  result  in  numerous  bonds  being  attached  to  the 
large  matrix  elements,  thereby  actually  destroying  the  pleasing  pattern.  An  example  of  this 
phenomenon  is  given  by  the  case 
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Note  that  row  and  column  permutations  can  transform  d*  into  d.  Since  d  is  in  block  form,  it 
conveys  more  information  about  the  group  structure  and  is  preferable  to  d*.  However,  if  any 

l.  TTili  procedure  can  he  modified,  however,  if  the  geac/tliasd  M«£  include*  only  tow-bond*  and  column-bond*  (aft  - 
not  necessarily  neared  neighbor),  and  lacks  diagonal  bonds. 
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of  the  three  following  MEs  (which  allow  next-to-nearest  neighbor  diagonal  bonds)  are  used, 
then  d'  is  preferable  to  d.  Evidently  optimization  of  the  bond  energies  leads  to  the  attachment 
of  as  many  (diagonal)  bonds  to  the  “500”  as  possible,  even  at  the  expense  of  block  form.  By 
contrast,  when  the  diagonal  bonds  are  excluded,  optimization  of  the  ME  will  produce  block 
form  if  this  is  possible.  (See  Appendix  B.) 


The  three  MEs  are 

ME,  (b)  =  £  6y  ay 
ij 

ME2  (b)  =  £  ay  ay 

« 

ME3(b)=£ay  0y 


where  5y  = 


1  by  >  0 

0  by  =  0 


aij =  “ij  +  l  +  ai,M +  ai+l,j  +ai-l,j 

+  l/2[a;+1J  +  1  +ai+1J.,  +ai.1J+1  ♦ai.,J.1] 

=  6i,j+l  +6i,j-l  +  5i+l,j  +  6i-  1.  j 

+  1/2  fo+l,j  +  l  +  5i+l,i-  1  +  5i- l.J  +  J  +  6  i- l.j-  J 

If  any  of  the  three  are  used,  then  d*  has  a  higher  ME  than  d: 

ME,  (d)  -  510  ME,  (d1)  =  1002 

ME2(d)=510  ME2  (d1)  =  2000 

ME3(d)  =  260.5  MEJ(d1)=1002 


C.  ADDITIONAL  PROPERTIES  OF  THE  MEs 


(l)  All  three  algorithms  are  unaffected  if  all  the  matrix  elements  are 
multiplied  by  a  positive  constant. 


(2)  /dl  three  algorithms  are  affected  if  a  constant  is  added  to  the  matrix 
dement.  In  particular  this  implies  sensitivity  to  the  choice  of  origin  of 
the  ordinal  scale  (e.g.,  0,1,2  versus  1,2,3)  when  rankings  are  used  as  the 
matrix  elements. 

(3)  The  sensitivity  of  the  Bond  Energy  Algorithm  to  the  choice  of  k 
depends  on  the  relative  magnitudes  of  the  various  matrix  elements.  If 
all  the  matrix  elements  are  0  or  1 ,  then  the  ME  is  independent  of  the 
choice  of  k.  If,  however,  the  matrix  elements  vary  greatly  in  magnitude, 
it  is  recommended  that  k  be  set  equal  to  2  instead  of  1.  This  choice 
preserves  scale  by  not  overemphasizing  the  largest  elements.  For 
example,  with  k  =  2,  the  bond  strength  between  elements  of  sizes  5  and 
7  will  be  the  \f  35,  close  to  their  average,  rather  than  the  inflated  value 
35  when  k  =  1. 
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A  MEASURE  OF  EFFECTIVENESS 
FOR  THE  MOMENT  ORDERING  ALGORITHM. 


It  has  been  pointed  out  that  one  of  the  properties  of  the  algorithm  is  to  drive  the  array 
being  operated  upon  toward  a  more  diagonal  form.  This  property  has  been  utilized  to  define1 
a  correlation  coefficient,  R,  to  measure  the  progress  of  the  iterative  procedure  and  the  quality 

of  the  final  result.  The  coefficient  has  been  defined  as  follows: 
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=  element  in  i^1  column  and  j1*1  row. 


t.  SugfettioQ  due  to  Dt.  Phillip  Gouid  of  IDA. 


Preceding  page  blank 


141 


Note  that  R  is  normalized  so  that  its  value  always  lies  between  zero  and  one.  For  the  special 
case  of  a  square  matrix  R=1  corresponds  to  only  the  main  diagonal  being  filled,  R=0  to  a 
random  distribution  of  value  throughout  the  array,  and  R=  -1  to  the  opposite  diagonal  only 
being  filled.  Initial  values  of  R  for  arrays  therefore  are  generally  near  zero,  and  as  the  algorithm 
proceeds  toward  a  solution,  R  generally  increases.  The  final  value  of  R  is  a  measure  of  the 
degree  of  diagonality  obtained  by  the  algorithm.  It  should  be  noted  however,  that  the 
algorithm  is  not  a  direct  attempt  to  maximize  R,  and  that  there  are  occasional  cases  in  which 
an  iteration  of  the  algorithm  will  decrease  R  instead  of  increasing  it. 
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APPENDIX  I 


MULTIPLE  SOLUTIONS  TO  THE  MOMENT  ORDERING  ALGORITHM 
FOR  A  SAMPLE  3x3  ARRAY 


MULTIPLE  SOLUTIONS  TO  THE  MOMENT  ORDERING  ALGORITHM 
FOR  A  SAMPLE  3x3  ARRAY. 


In  order  to  investigate  the  factors  which  lead  to  multiple  solutions  to  the  Moment 
Ordering  Algorithm,  the  following  experiment  was  carried  out.  It  deals  with  a  3  x  3  array,  hut 
it  is  believed  that  the  conclusions  drawn  may  be  useful  in  understanding  the  phenomenon  for 
the  vastly  more  complicated  cases  of  larger  arrays. 

1.  A  sample  3x3  array  (Fig.  1*1)  was  constructed.  For  simplicity,  its  rows  were 
each  nomalized  to  10.  Two  of  the  rows  were  Fixed  (7.2.1  and  3,5,2),  while  the 
elements  in  the  third  were  allowed  to  take  on  various  valu»  ..  (always  subject  to 
the  normalization  and  the  restriction  that  all  elements  be  non-negative). 
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FIGURE  1-1.  Experimental  3x3  Array 

2.  For  every  possible  combination  of  values  for  the  elements  of  the  third  row.  the 
resulting  array  was  analyzed.  In  particular,  the  number  of  possible  stable 
solutions  was  determined. 

3.  The  results  are  presented  in  Fig.  1-2.  Every  point  inside  the  triangle  represents  a 
possible  third  row  of  the  array.  The  values  of  the  three  elements  are  read 
upwards  from  each  face.  (Note  that  the  sum  of  the  distance  from  any  point 
inside  to  all  three  faces  of  the  equilateral  triangle  is  constant -in  this  case  equal 
to  10.)  The  sets  of  elements  corresponding  to  the  First  two  rows  are  marked  as 
A  and  B,  and  for  each  other  point  the  multiplicity  of  solutions  to  the  resulting 
array  is  shown. 

4.  The  resulting  overall  pattern  indicates  that  when  the  third  point  iscolinear.  or 
nearly  so,  with  points  A  and  B  (that  is,  when  the  three  rows  are  in  a 
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well-denned  order  within  the  coordinate  system  ol  rig.  1-2)  only  or 
usually  exists.  However,  in  the  region  in  the  lower  left  hand  secti 
triangle,  where  the  third  point  forms  a  triangle  with  A  and  B,  rat 
straight  line,  three  stable  solutions  exist.  This  indicates  that,  when  o 
linear  ordering  exists,  the  algorithm  will  find  that  ordering,  but 
several  orderings  are  equally  satisfactory,  it  may  find  each. 


ONE  SOLUTION 

f . "1  TWO  SOLUTION 

H  THREE  SOLUTIOl 


FIGURE  1-2.  Multiplicity  of  Solutions  for  a  Small  Array 
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